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INTRODUCTION 


Background 

The  ultimate  objective  of  this  effort  is  to  study  the  hypersonic  flow 
field  associated  with  a  shaped  charge  jet.  Of  concern  will  be  how  and  to 
what  extent  perturbations  of  the  jet  are  caused  by  interactions  of  the  flow 
field  with  various  target  geometries.  This  report  details  the  first  part  of 
the  effort  to  develop  an  axisymmetric  fluid  dynamics  code:  the  derivation  of 
applicable  equations  and  the  development  of  a  one  dimensional  code.  Many 
basic  relationships  applicable  to  both  codes  are  derived  here. 

Considerable  research  has  been  conducted  in  the  fields  of  jet  formation 
and  jet  penetration,  but  little  effort  has  been  devoted  to  studying  the 
aerodynamic  forces  that  influence  the  jet.  At  typical  jet  velocities  of 
greater  than  Mach  20,  an  extremely  strong  flow  field  is  developed  which  may 
cause  severe  perturbations  in  the  jet.  For  example,  an  extremely  strong  bow 
shock  wave  is  formed.  As  the  jet  flies  through  the  hole  in  a  guidance 
package,  or  through  the  hole  in  a  skirting  plate  made  by  the  leading  portion 
of  the  jet,  the  bow  shock  is  reflected  back  to  impinge  on  a  trailing  portion 
of  the  jet,  applying  potentially  disruptive  forces.  Another  potentially 
disruptive  force  results  from  the  strong  wake  behind  a  jet  particle.  The 
subsequent  particles  must  fly  through  this  turbulence.  The  situation  may  be 
even  more  severe  for  particles  farther  back  in  the  jet. 

The  problem  is  studied  by  numerical  solution  of  the  equations  of  motion 
utilizing  a  method  developed  by  S .  K.  Godunov,  a  Russian  mathematician. 

This  first  order  inviscid  method  is  suitable  to  hypersonic  flow  regimes.  It 
has  been  used  by  the  Launch  and  Flight  Division  of  the  BRL  to  solve  a 
geometrically  similar  problem  of  a  projectile's  flight  through  a  muzzle 
brake.  This  problem,  however,  involved  only  one  projectile  at  less  than 
hypersonic  velocities.  The  code  will  be  validated  by  solving  some  problems 
with  known  solutions.  Then  a  cylindrical  hole  of  circular  cross  section 
will  be  studied.  These  calculations  will  be  experimentally  verified. 
Finally,  more  complex  geometries  may  be  numerically  studied. 

In  the  first  phase  of  .his  study,  reported  here,  the  necessary  equations 
are  derived  and  a  one -dimensional  code  is  written  for  an  inviscid,  perfect 
gas.  Theoretical  and  experimental  verifications  are  carried  out  which  show 
that  the  technique  accurately  predicts  those  flow  phenomena  critical  to  this 
study:  shock  and  expansion  wave  formation,  wave  reflection,  and  wave /wave 
collisions . 


The  Godunov  Technique 

The  one -dimensional  scheme  developed  by  Godunov  is  explained  and 
expanded  to  two  dimensions  (and  ax i- symmetric)  in  a  1961  paper  in  the 
Russian  Journal  of  Mathematics  and  Mathematical  Physics  (Ref  1) .  The 
principal  advantage  of  the  method  lies  in  the  physically  realistic  approach 
used  to  obtain  the  mathematical  solution  of  the  equations  of  fluid  motion. 
That  the  solution  technique  can  be  so  well  expressed  in  physical  terms  is 
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fascinating,  and  a  fine  comment  on  the  beauty  with  which  mathematics 
expresses  the  physical  laws  of  nature.  In  fact,  Godunov  (Ref  2)  derives  all 
the  necessary  equations  (except  the  conservation  equations)  purely  from 
mathematical  considerations.  As  a  final  step,  he  indicates  that  his 
equations  are  the  same  as  those  that  would  result  from  consideration  of  a 
physical  system.  The  physical  system  will  be  the  basis  of  the  equations 
derived  here. 

The  physical  system  referred  to  above  is  the  one  dimensional  fluid 
discontinuity,  sometimes  called  the  Riemann  problem  (after  the  German 
mathematician  who  was  the  first  to  attempt  to  calculate  shock  properties)  or 
the  shock  tube  analogy.  To  see  how  this  Riemann  problem  helps  in  the 
solution  of  the  equations  of  fluid  motion,  consider  a  one  dimensional  fluid, 
divided  into  finite  ceils  and  at  initial  time  t0  containing  the  pressure 
distribution  as  shown  in  figure  la.  Similar  curves  can  be  imagined  for  the 
other  dependent  variables,  velocity  (u) ,  density  (p),  internal  energy  (e). 
The  problem  being  addressed  here  is  one  of  unsteady  flow,  so  these 
quantities  will  be  functions  of  time  as  well  as  functions  of  x. 

Thus  the  gas  j.s  divided  into  one-dimensional  regions  of  thickness  Ax. 

Some  appropriate  average  values  for  p,  u,  p,  e  are  chosen  for  each  layer 

(at  say,  t  -  r).  These  are  denoted  by  p  ,  u  ,  etc.  Consequently, 

m- 1/2  m-i/2 

the  values  in  neighboring  layers  may  not  be  the  same,  as  seen  in  figure  lb. 
Because  of  the  incremental  nature  of  the  assumed  form  of  the  pressure 
distribution,  a  pressure  discontinuity  exists  at  the  boundary  between  cells. 
This  is  analogous  to  the  classical  shock  tube  problem,  which  considers  the 
solution  of  the  Riemann  problem,  as  depicted  in  figure  lc.  Of  course  in  the 
present  case  the  diaphragm  is  imaginary.  It  is  imagined  to  rupture  by 
unfreezing  the  time  variable  and  permitting  the  discontinuities  in  pressure, 
density  and  velocity  to  seek  equilibrium.  The  discontinuities  at  x^  are 

thus  resolved  into  a  shock  wave,  expansion  wave  and  contact  discontinuity 
(or  some  other  combination) ,  propagating  from  the  cell  boundary  at  x^  ,  or 

simply  m,  toward  its  neighbors  at  m-1  and  m+1. 

The  conditions  behind  these  waves  are  well  known  from  shock  tube  theory, 
and  can  be  computed  from  the  known  initial  conditions  on  the  two  sides  of 
the  diaphragm.  More  importantly,  the  conditions  behind  the  waves  are 
constant.  The  essence  of  the  Godunov  technique  is  to  utilize  these  constant 
conditions  at  the  cell  boundaries  to  compute  the  flux  of  properties  into 
each  cell  during  the  time  step.  The  conservation  equations  are  used  to 
compute  new  average  properties  in  each  cell  and  the  entire  process  is 
repeated  in  the  next  time  step.  The  integral  form  of  the  conservation 
equations  must  be  used  if  the  solution  is  to  be  valid  for  flow  fields 
containing  discontinuities  (shock  waves). 

The  equations  which  constitute  the  general  solution  of  the  Riemann 
(shock  tube)  problem  are  nonlinear.  For  weak  (sonic)  waves,  the  equations 
can  be  linearized.  Figures  la  and  lb  show  that  the  magnitude  of  the 
discontinuities  in  each  individual  Riemann  problem  (which  result  from  the 
approximation  to  the  continuous  distribution)  are  governed  in  part  by  the 
size  of  the  increment  Ax.  Thus  the  resulting  waves  can  be  guaranteed  to  be 
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weak  waves  by  selection  of  a  sufficiently  fine  grid.  The  consequential 
linearized  Riemann  system  is  valid  in  all  cases  except  the  one  in  which  a 
discontinuity  in  the  flow  field  (e.g.  bow  shock  wave,  reflected  wave,  etc.) 
actually  passes  through  one  of  the  two  cells. 

To  insure  validity  of  the  Riemann  model  when  approximating  continuous 
flow,  the  fluid  properties  at  the  point  x  must  remain  constant  during  the 

time  step;  they  must  not  be  disturbed  by  the  waves  propagating  from  the 
neighboring  points  x  and  x  .  That  is,  the  properties  at  x  must  only 

m- 1  tTT+,1  m 

result  from  the  resolution  of  the  discontinuity  at  x^.  The  fact  that  the 

region  between  points  x  and  x  (or  x  and  x  )  is  complicated  by  the 

m  m+1  m  m-1 

presence  of  additional  waves  emanating  from  x  (or  x  )  is  of  no 

fTHi  m-  1 

consequence  to  the  properties  at  x  .  The  properties  computed  for  point  x 

m  fn 

will  remain  valid  until  the  waves  from  the  neighboring  points  arrive  at  x  . 

m 

Thus  the  time  interval  At  must  be  less  than  the  time  required  for  the  wave 
to  traverse  che  distance  Ax.  If  W  is  the  wave  speed  relative  to  the  fixed 
coordinate  frame , 


At  < 


Ax 

W 


is  the  physically  required  relationship  between  the  step  sizes.  Godunov 
also  derives  this  criterion  from  a  purely  mathematical  stability  analysis, 
again  demonstrating  the  harmony  between  physical  and  mathematical 
descriptions  of  nature. 


General  Approach  to  Solving  the  Fluid  Dynamics  Problem 

The  Godunov  technique  requires  the  solution  of  two  fluid  dynamics 
problems.  The  primary  problem  involves  use  of  the  conservation  laws  to 
update  the  properties  in  each  cell  based  on  the  flux  of  each  property  across 
the  cell  boundaries.  Since  dicontinuities  may  exist  in  the  flowfield,  the 
integral  form  of  the  conservation  equations  must  be  used. 

The  secondary  fluid  dynamics  problem  to  be  solved  is  the  Riemann  problem 
at  the  cell  boundaries.  The  fluid  condition  behind  the  waves  constitutes 
the  solution,  and  these  conditions  are  used  in  the  conservation  laws  (the 
primary  problem)  to  determine  fluxes  and  hence  new  average  cell  properties. 

Figure  2  shows  a  schematic  of  the  conditions  in  the  Riemann  problem. 

For  simplicity,  the  subscript  k  denotes  conditions  in  cell  (m-1/2)  and 
subscript  1  denotes  conditions  in  cell  (m+1/2) .  The  negative  running  wave 
and  positive  running  wave  may  also  fall  in  the  same  quadrant  (either  one) , 
but  by  definition  the  positive  running  wave  has  the  more  positive  absolute 
(i.e.  relative  to  the  fixed  cell  boundary)  velocity.  The  location  of  each 
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wave  (left  or  right  quadrant)  will  depend  on  the  initial  fluid  velocities  in 
the  two  cells,  and  the  resultant  wave  velocity  relative  to  the  fluid. 

From  classical  shock  tube  theory,  the  pressure  and  velocity  of  the  gas 
on  either  side  of  the  contact  discontinuity  are  known  to  be  equal  and 
constant;  hence  the  terminology  p0  and  u0  in  the  regions  behind  the  waves, 
regions  2  and  3.  Similarly,  the  density  is  known  to  be  different  on  the  two 
sides  of  the  contact  discontinuity,  but  to  be  constant  in  each  of  the 
regions  2  and  3.  That  is,  the  density  is  discontinuous  across  the  contact 
discontinuity.  Pi  ,  uj  ,  px  ,  and  p4  ,  u4  ,  p4  are  the  known  and  constant 

conditions  of  the  gas  prior  to  the  passing  of  the  wave.  Note  that  a 

compressive  wave  is  characterized  by  a  higher  pressure  behind  the  wave  than 
in  front  (  p0  >  Pi  or  p0  >  p4  )  while  an  expansion  wave  leaves  a  lower 

pressure  behind  (  p0  <  Pi  or  Po  <  p4  ) • 

Since  the  integral  equations  must  be  used  for  discontinuties  in  the  flow 
field,  they  will  be  used  to  derive  the  conditions  behind  the  shock  wave.  The 
expansion  wave  is  not  a  discontinuity  in  the  flow,  so  the  less  awkward 
differential  form  of  the  equations  of  fluid  motion  will  be  used  to  derive 
the  conditions  behind  the  expansion  wave. 
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Chapter  1 


THE  EQUATIONS  OF  1NVISCID  FLUID  MOTION 


The  conservation  equations  (or  equations  of  motion)  are. used  three  times 
in  the  following  chapters:  in  Chapter  2  for  the  derivation  of  the  finite 
difference  equations  that  are  used  to  compute  the  flux  of  properties  into 
the  cell  during  the  time  step;  in  Chapter  3  to  derive  the  shock  wave 
equations  for  the  Riemann  problem;  and  again  in  Chapter  3  to  derive  the 
expansion  wave  equations  for  the  Riemann  problem. 

To  derive  these  conservation  laws  using  a  Lagrangian  approach  is 
aesthetically  more  satisfying  than  using  an  Eulerian  approach.  In  the 
Lagrangian  approach,  a  particular  element  of  fluid  is  studied  as  it  flows. 
The  control  volume  in  this  case  is  that  which  encloses  the  element  of  fluid 
under  study.  It  can  change  size  and  shape  as  it  moves  along  in  the  flow 
field,  but  unlike  the  Eulerian  control  volume,  the  mass  within  it  is 
constant.  For  this  reason  the  Lagrangian  approach  seems  more  naturally 
suited  to  the  derivation  of  conservation  laws. 

Both  the  differential  and  the  integral  form  of  the  conservation 
equations  are  derived  in  this  chapter.  To  arrive  at  the  differential  form 
of  the  conservation  equations,  total  derivatives  of  volume  integrals  are 
expressed  in  Eulerian  terms  (the  differentiation  is  moved  inside  the 
integral)  by  use  of  the  Reynold's  Transport  Theorem.  This  form  is  then  used 
to  generate  the  finite  difference  equations  (which  compute  the  flux  of 
properties  into  each  cell  during  the  rime  step),  and  to  derive  the  expansion 
wave  equations  for  the  Riemann  problem.  The  derivation  of  the  shock  wave 
equations  for  the  Riemann  problem  is  facilitated  by  use  of  the  integral  form 
of  the  conservation  equations.  This  form  results  directly  from  a  Lagrangian 
derivation,  as  detailed  in  the  latter  portion  of  this  chapter. 

The  Differential  Equations  of  Motion 

The  Reynold's  Transport  Theorem  is  required  for  the  derivation  of  the 
differential  equations  of  motion.  In  deriving  the  Reynold's  Transport 
Theorem,  use  is  made  of  the  Gauss  Divergence  Theorem.  It  relates  a  surface 
integral  to  a  volume  integral  and  is  written  (ref  3) : 

h*wdcr  —  V*wdr 

Js  Jv 


where  w  is  any  vector  function  of  time  and  position,  da  is  a  scalar 
element  of  surface  area,  h  is  a  unit  vector  normal  to  da  and  directed 
outward,  and  dr  is  a  scalar  element  of  volume.  Replacing  hda  with  dS 
and  dr  with  dV  (also  a  scalar) ,  the  Gauss  Divergence  Theorem  is  written 

I  w*dS  «•  I  V*cl>  dV 

J  e  J  w 
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For  a  scalar  function  of  time  and  position  (w) ,  the  Gauss  Divergence  Theorem 
is  (ref  3,  eq.  9.124) 


I  u>  dS  -  Vu  dV 

Js  Jv 


(102) 


Reynold's  Transport  Theorem 

In  the  Lagrangian  approach  (figure  3),  the  moving  mass  of  fluid  being 
studied  is  imagined  to  be  enclosed  by  a  moving  control  volume  that  may 
change  size  and  shape,  but  always  contains  the  same  mass;  i.e.  no  mass 
cresses  its  surface.  An  observer  moving  with  this  Lagrangian  control  volume 
can  detect  no  change  in  a  fluid  property  a  with  respect  to  any  spatial 
coordinates,  since  to  him  the  control  volume  and  any  mass  within  it  are 
stationary.  But  the  property  a  is  changing  as  time  passes,  so  in  the 
Lagrangian  sense,  a  -  a( t)  only. 

Of  interest  is  the  integral  of  the  property  a  over  the  volume  V: 


I(t)  -  I  a(t)  dV 
•'V(t) 


Furthermore,  the  time  rate  of  change  of  this  integral  is  desired: 


dlft) 

dt 


dt  f 


o(t)  dV 
V(t) 


(103) 


(104) 


The  fundamental  theorem  of  the  differential  calculus  is  used  to  expand 
eq.  104  to  the  following  form. 

k  I  ■><t)  dV  -  its  k  I  f  dV  - 

JV(t)  L  JV<t+4t) 

The  geometry  of  figure  3  is  used  to  modify  this  expression  to  the  following 
intermediate  form  of  the  Reynold's  Transport  Theorem. 


1 


a(t) 

V(t> 


dV 


(105) 


d_ 

dt 


a(t)  dV  - 


da(t) 

dt 


dV  + 


I  a(t)  fvdS] 
J  set) 


(106) 


Note  that  the  time  derivative  of  the  quantity  I(t)  (eq.  103)  at  a  certain 
time  t  has  been  expressed  in  terms  of  quantities  evaluated  at  the  same  time. 
This  has  important  implications  for  the  explicit  finite  difference  technique 
being  derived. 
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The  Gauss  Divergence  Theorem  (eq.  101)  converts  the  surface  integral 
into  a  volume  integral,  and  eq.  106  becomes  the  following. 


dt  I 


a( t)  dV 


f  r  2s 

J  at 

JV(t)  L 


+  V*QV 


(107) 


Equation  107  is  the  Reynold's  Transport  Theorem.  It  accomplishes  two 
things.  Mathematically,  it  moves  the  differentiation  inside  the  integral. 
Physically,  it  relates  Lagrangian  terms  on  the  left  hand  side  to  Eulerian 
terms  on  the  right  hand  side. 


Conservation  of  Mass 

In  a  Lagrangian  frame  of  reference,  the  mass  within  the  control  volume 
does  not  change.  An  element  of  mass  within  V  is  written  p dV.  The  total 
mass  of  fluid  within  V  is 

m  —  T  pdV 
•*  vet) 


Since  the  mass  does  not  change  with  time,  the  time  derivative  of  the 
integral  must  be  zero. 


dm 

dt 


pdV  -  0 


(111) 


With  the  Reynold's  Transport  Theorem  (eq.  107),  eq.  Ill  becomes 

ttf  <’<t)dV-f  [ff  *  V-p^ldV  -  0 
JV(0  JV<t)  L  J 

and  since  the  volume  is  arbitrary  (see  also  ref  4,  pg  153),  the  only  way  for 
the  integral  to  be  zero  for  all  control  volumes  V  is  if  the  integrand  itself 
is  zero  at  every  point  in  V. 


y.pv 


0 


(112) 


Equation  112  is  the  general  differential  form  of  the  continuity  equation  in 
the  three  dimensional  case.  For  a  one  dimensional  problem,  it  reduces  to 


(p  u) 


0 
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(113) 


Conservation  o£  Momentum 

The  resultant  force  acting  on  the  fluid  within  the  control  volume  is 
equal  to  the  time  rate  of  change  of  momentum. 


£F  - 


dM 

dt 


(114) 


Neglecting  body  forces  (gravitational,  electro-magnetic)  and  viscous 
shear  forces  (inviscid  fluid) ,  the  only  force  acting  on  the  fluid  within  V 
is  due  to  the  pressure  acting  on  the  surface  of  V. 


p(-dS) 


(115) 


Here  the  negative  sign  is  required  since  dS  is  positive  outward,  but  the 
resultant  force  is  to  be  positive  inward. 


The  fluid  within  the  control  volume  V  changes  momentum  due  to  changes  in 
velocity  resulting  from  its  unsteady  motion  or  its  flowing  around  objects. 
The  momentum  of  an  element  of  this  fluid  is  (pdV)v.  The  momentum  of  the 
entire  volume  of  fluid  is 


M 


I 


pv  dV 


Substituting  this  and  eq.  115  into  eq.  114, 


-  f  P  dS  -  ^  f  pv  dV 
J  c  J  \/ 


The  surface  integral  may  be  converted  to  a  volume  integral  by  use  of  the 
scalar  form  of  the  Gauss  Divergence  Theorem,  eq.  102. 

-  I  Vp  dV  -  Jj-  f  pv  dV 
J  u  Ju 

(116 

The  integral  on  the  right  hand  side  of  eq.  116  can  be  converted  to  Euler ian 
terms  by  the  Reynold's  Transport  Theorem  (eq.  107).  Since  a(t)  in  eq.  107 
is  a  scalar,  the  quantity  pv  in  eq .  116  must  first  be  written  in  scalar 
form.  That  is,  the  vector  momentum  equation  will  have  to  be  rewritten  in 
its  three  equivalent  scalar  components  first.  For  the  present  work,  only 
the  x-momentum  equation  is  of  interest.  Thus  eq .  116  becomes 
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-  1  SdV  -  dH  I  dv 

With  the  Reynold's  Transport  Theorem  (eq.  107),  this  becomes 

J  [  ^(pu>^f-V*puv  +  j  dV  -  0 


Since  V  is  arbitrary, 


+  V • puv  + 


d£ 

3x 


0 


(117) 


Equation  117  is  the  general  differential  form  of  the  x-momentum 
equation.  Note  that  it  is  a  scalar  equation.  For  a  one  dimensional 
problem,  it  reduces  to 


k  (pu)  +  k  (pu2) 


& 

3x 


(118) 


This  can  be  expanded  to  the  following  form. 

.  iL  ,...s  1  +  &u  £u  l3£ 

at  +  3x  ^pu)  J  +  3t  +  U  3x  p  3x 


ij 

P 


But  the  term  in  brackets  is  zero  by  the  continuity  equation  (eq.  113). 
Thus  the  one  dimensional  momentum  equation  is: 


3u  3u 

+  u  T“ 
3t  3x 


1 

p  3x 


(119) 


Conservation  of  Energy 

The  principle  of  conservation  of  energy  states  that  the  rate  of  change 
of  energy  of  the  fluid  within  the  control  volume  V  is  equal  to  the  rate  at 
which  heat  is  added  plus  the  rate  at  which  work  is  done  on  the  fluid.  The 
energy  balance  is  written: 


rate  at  which  heat 
is  added  to  fluid 
within  V 


rate  at  which  work 
is  done  on  fluid 
within  V 


rate  of  change  of 
energy  of  fluid 
within  V 


(120) 


Heat  conduction  and  radiation  in  the  fluid  will  be  neglected.  This  is 
permissible  for  the  very  short  times  being  computed  (milliseconds) .  Since 
no  mass  enters  or  leaves  the  Lagrangian  control  volume,  no  heat  can  be 
convected  across  S,  the  surface  of  V.  Thus  the  first  term  in  the  energy 
balance  is  zero. 


9 


The  work  in  the  second  terra  is  done  by  forces;  since  body  forces  and 
shear  stresses  are  being  neglected,  the  only  active  force  is  due  to  the 
pressure  on  S,  pdS.^  The  rate  of  work  done  by  this  force  on  a  fluid  moving 
at  velocity  v  is  pdS*v.  The  total  rate  of  work  done  by  the  pressure 
exerted  on  the  fluid  within  V  is 


W-.  f 


pv«  dS 


(121) 


where  the  negative  sign  is  required  to  make  W  positive  for  pressure  acting 
inward,  since  p  and  dS  are  positive  in  opposite  directions. 

The  specific  energy  of  a  mass  of  fluid  is  its  specific  internal  energy 
e,  plus  its  specific  kinetic  energy,  v2/2,  where  v2  -  v*v,  a  scalar. 

All  other  forms  of  energy  (potential,  chemical,  etc.)  have  been  neglected. 

Since  the  mass  of  an  element  of  fluid  within  V  is  pdV  and  its  specific 
energy  is  (e  -f  v2/2),  the  total  energy  of  the  fluid  within  V  is 


E  -  [  p[e  +  v2/2 ]  dV 

J  V 


(122) 


Substitute  eq's.  121  and  122  into  eq.  120,  recalling  that  the  first  term  in 
eq.  120  is  zero. 


-  f  pv*dS  -  ~  j  p[e  +  v2/2]  dV 

J  q  J  u 


The  term  on  the  left  hand  side  is  converted  to  a  volume  integral  by  the 
Gauss  divergence  Theorem  (eq.  101). 


-  J  V«pv  dV 

Jv 


d_  r 
L 


p [ e  +  v2/2 ]  dV 


The  term  on  the  right  hand  side  is  simplified  by  the  Reynold's  Transport 
Theorem  (eq.  107).  After  rearranging,  the  following  equation  results. 


J.[fc 


p[e  +  v2/2]  +  V*p[e  +  v2/2]v  +  V*pv  dV  -  0 


Since  V  is  arbitrary, 


p[e  +  v2/2]  +  V«p[e  +  v2/2]v  +  V*pv 


(123) 
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Equation  123  is  the  general  differential  form  of  the  conservation  of  energy. 
For  a  one  dimensional  problem,  it  reduces  to  the  following. 


4~  p[e  +  u2/2]  +  ~  pu[e  +  u2/2  +  ^]  =  0 

P  (124) 


The  Integral  Equations  of  Motion 

In  the  Lagrangian  approach,  we  consider  a  certain  body  of  fluid  as  it 
moves  along  in  the  flowfield.  The  integral  equations  are  simplified  by  this 
approach  since  the  mass  of  the  fluid  element  under  study  does  not  change 
with  time.  As  seen  in  figure  4,  this  fluid  occupies  a  region  of  the  x  axis 
from  a0(t)  to  aj(t),  and  contains  a  discontinuity  (shock  wave)  at  x  =  £(t). 
Note  that  subscript  0  denotes  conditions  behind  the  shock  wave,  and 
subscript  1  denotes  conditions  ahead.  This  is  consistent  with  figure  2. 


Conservation  of  Mass 

In  this  one  dimensional  analysis  we  consider  a  unit  area  in  the 
direction  of  flow  (i.e.  perpendicular  to  the  x  axis).  The  mass  of  an 
element  of  this  fluid  is  then 

dm  -  p(x, t )  [ldx] 

and  the  total  mass  of  this  body  of  fluid  is 

i<” 

m  -  p  ( x ,  t )  dx 

Ja0(t) 


The  mass  of  this  fluid  is  constant;  this  is  implicit  in  the  Lagrangian 
approach.  The  principle  of  conservation  of  mass  follows  directly: 


ax(t) 
a  0?  t  > 


p(x,t)  dx  -  0 


(125) 


Conservation  of  Momentum 

The  momentum  of  an  element  of  the  fluid  is 
dM  -  p(x,t)  [ldx]  u(x,t) 

so  the  total  momentum  of  the  body  of  fluid  is: 
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pai(t) 

M  -  p<x,t)  u(x,t)  dx 

Ja0(t) 


The  only  forces  acting  on  our  one-dimensional  fluid  mass  are  the 
pressures  at  the  ends,  l*p0  and  -l«p1(  where  1  is  the  area  of  each  end. 

Body  and  shear  forces  are  neglected.  Combining  the  fluid  momentum,  pressure 
forces,  and  Newton's  Law  yields  the  principle  of  conservation  of  momentum: 


“  f 

dt  J 


■a  O 
SqCO 


p<x,t)  u(x,  t)  dx  =  p0 


(126) 


Note  that  both  velocity  and  pressure-force  have  been  considered  positive  in 
the  +x  direction. 


Conservation  of  Energy 

If  the  body  of  fluid  is  not  undergoing  any  chemical  reactions  and  has  no 
heat  transferred  in  or  out,  we  may  assume  it  changes  energy  only  because  of 
the  work  done  by  the  external  forces,  l*p0  and  l*Pi .  The  rate  at  which  work 
is  done  (power)  must  be  equal  to  the  rate  at  which  the  energy  of  the  body  of 
fluid  is  increasing.  This  is  the  conservation  of  energy. 

p0  does  positive  work  at  the  rate  of  l*p0u0 
px  does  negative  work  at  the  rate  of  l*p1u1 

If  e(x,t)  is  the  specific  internal  energy  and  u(x,t)2/2  is  the  specific 
kinetic  energy,  and  all  other  forms  of  energy  are  neglected  (chemical, 
potential,  etc.),  then  the  total  energy  of  an  element  of  the  fluid  is 


dE  -  [e(x,t)  +  u(x,t)2/2]  p(x,t)  [ldx] 


and  the  total  energy  of  the  entire  fluid  mass  is 
ra1(t) 

E  -  (e(x,t)  +  u(x,t)2/2 ]  p<x,t)  [ldx] 

Ja0<'> 


Finally,  the  conservation  of  energy  follows. 


d_ 

dt 


ra^t) 

[e(x,t)  +  u(x,t)2/2] 
Ja0(t) 


p(x,t)  dx 


Pouo 


PlU! 


(127) 
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Chapter  2 


THE  FINITE  DIFFERENCE  EQUATIONS 


The  following  differential  equations  of  motion  for  one -dimensional 
inviscid  fluid  flow  problems  were  derived  previously. 


d_ 

at 


3_ 

at 


[p] 

+  L 

[pu] 

-  0 

[pu] 

+  k 

[p+pu2] 

-  0 

a_ 

at 


u 


P(e  +  *-) 


3x 


,  U2  D . 

P  u(e  +  —  + 


(113) 

(118) 

(124) 


These  equations  are  valid  at  any  point  in  the  flow  field  (and  at  any  time) 
provided  the  flow  variables  are  changing  continuously  at  that  point  (and 
time).  If  points  (times)  exist  in  the  flowfield  at  which  the  flow  variables 
are  changing  discontinuously ,  the  differential  equations  of  motion  must  be 
integrated  over  some  arbitrary  but  finite  area  containing  the  discontinuity. 
In  this  way  the  discontinuity  is  made  mathematically  continuous. 

For  example,  consider  the  pressure  discontinuity  in  space  (i.e.  shock 
wave)  shown  in  figure  5.  The  derivative  dp  is  clearly  infinite  at  x=b , 

owing  to  the  pressure  discontinuity  at  that  point.  This  difficulty  would 
made  an  equation  like  eq.  118  meaningless.  However,  if  the  derivative  is 
integrated  over  some  finite  interval  containing  the  discontinuity,  the 
result  is  well  behaved.  The  discontinuity  can  be  written  mathematically 
with  the  help  of  the  Heavyside  step  function  (ref  3) . 


H(x-b) 


x  >  b 
x  <  b 


Thus  the  pressure  function  pictured  in  figure  5  is  written  as  follows. 


p(x)  -  Pi  +  [p2  -  Pi  ]  H(x-b) 


The  derivative  of  this  discontinuous  function  is  of  concern.  The  derivative 
is  integrated  over  a  region  containing  the  discontinuity,  say  x=a  to  x=c . 
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r  s dx  ■  j°  dx  p1+(p2-pi)H<)<-b>]dx 


(p 


:-p>)  f  H 


'(x-b)  dx 


Notice  that  and  p2  are  constant  in  the  region  of  integration.  Since  the 
integral  of  a  derivative  of  a  function  is  simply  that  function,  the  integral 
above  is  simply  the  Heavyside  function  evaluated  at  x=c  minus  the  function 
at  x-a  (1  and  0  respectively,  by  definition) .  The  original  integral  is 
rewritten,  showing  that  the  discontinuity  has  been  mathematically  avoided: 


r 


dx  -  (p2-Pl) 


Since  discontinuities  in  the  flow  field  can  exist  in  either  time  or 
space,  the  differential  equations  of  motion  must  be  integrated  twice. 
Consider  an  arbitrary  flow  variable  ^(x,t)  (density,  momentum,  or  energy)  as 
shown  in  figure  6.  Notice  in  the  figure  that  the  letters  f  and  F  are  used 
to  indicate  the  value  of  the  fluid  property  <j>  on  the  boundaries  of  the  cell. 
These  will  be  discussed  in  greater  detail  below.  The  derivatives  of  the 
property  4>  are  integrated  over  the  x-t  region  shown  in  figure  6. 


n 


r-’  f  d,  dt 


f  3^ 

LJX  a* 

m 


dx 


dt 


<£(x  ,t) 

m 


dt 


(201) 


Here  <£(x  ,t)  is  the  instantaneous  value  of  <f>  at  the  point  x  ,  so  <f>(x  ,t)  is  a 
in  mm 

function  of  time  only.  It  is  now  assumed  to  be  constant  over  the  time  step 
(from  t  to  t  ),  although  it  may  vary  from  one  time  step  to  the  next.  A 

similar  assumption  is  made  about  ^(x  ,t).  Denote  these  values  bv  F  and  F 

m* 1  m  m+ 1 

These  two  values  are  defined  by  the  following  expressions. 


F 

m 


4>(x  ,t)  dt 
m 


(202) 


F 

m+1 


0(x  _,t)  dt 
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(203) 


That  is  F  is  a  time  averaged  value  (over  a  single  time  step)  taken  along  the 
beginning  or  ending  cell  boundary.  Henceforth,  all  upper  case  variable 
names  (R  for  density,  U  for  velocity,  E  for  energy)  will  represent  such  time 
averaged  values.  They  are  also  the  fluid  properties  behind  the  waves  that 
are  imagined  to  emanate  from  the  cell  boundaries  due  to  the  discontinuities 
in  fluid  properties  that  exist  there.  They  are  found  by  solving  the  Riemann 
problem  at  the  boundary,  which  is  detailed  in  chapter  3.  Substitution  of 
eqs.  202  and  203  into  eq.  201  yields  the  desired  result. 


rVi  rx 

J  t  "  X 


m+1 


dx  dt  -  r  F  .  -  F 

m 


(204) 


The  time  derivative  of  4>  is  integrated  in  an  analogous  manner. 


fXn>*1  l^n+l  dd)  ,  ,  rXm+1  f"  fCn+1  3^ 

Jx  Jc  LJC 


px  r  i 

^  b(x,t  )  -  <j>(x  ,t)  dt 

Jx  L  n  J 

m  (205) 

Here  ^(x.t^)  represents  the  value  of  <j>  at  any  point  at  the  time  t  ;  that 

is  t  )  is  a  function  of  x  only.  It  is  now  assumed  to  be  constant  over 

n 

the  space  step  (the  1-D  analogy  of  the  2-D  computational  cell)  x  to  x 

m  m+1 

although  <f>  may  vary  from  one  position  (space  step)  to  the  next.  A  similar 

assumption  is  made  about  <t>{x, t  ).  These  average  values  are  denoted  bv 

n+1 

f  and  f  (see  figure  6).  In  this  notation,  the  m+i/2  denotes  the 
m+1/2 

average  value  over  the  x  interval  from  m  (i.e.  x)  tom+i  (i.e.  x  ),  the 

m  m+1 

superscript  denotes  an  average  taken  at  the  end  of  the  time  step,  t  ,  and 
the  subscript  denotes  an  average  taken  at  the  beginning  of  the  time  step,  t  . 
The  average  value  of  a  property  may  be  written  in  the  following  manner. 


-  rx 

i  m+1  .  .  »  . 

*hj  dX 

J  v 


(206) 


(  207) 
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Thus  f  is  a  space  averaged  value  (over  a  single  space  step)  taken  at  either 
the  beginning  or  end  of  the  time  step.  Henceforth,  lower  case  variable 
names  ( p  for  density,  u  for  velocity,  and  e  for  energy)  will  be  used  to 
denote  the  fluid  properties  that  exist  within  the  cell.  They  have  been 
assumed  to  be  constant  within  the  cell.  Notice  that  the  spatial  averaging 
procedure  forces  any  discontinuities  in  the  flow  field  to  be  displaced  to 
the  border  of  a  computational  cell.  Substitution  of  eq .  206  and  eq.  207 
into  eq.  205  yields  the  desired  result. 


fXm+1  fCn+1  3  .  .  ,  fSi+l  f' 

-  [pu]  dt  dx  +  I 

Jx  J  t  J  t  J> 


t  -X 


m+1  5 


[p+pu2]  dx  dc  =  0 


“  x 

n  m 


Again  eqs .  204  and  208  are  used  to  reduce  this  to  the  finite  difference 
form. 


[pu]  „  -  ~  [P+RU2]  -  [P+RU2 

J m+l/2  h  L  J n>+i  L 


(210) 


Equation  210  relates  the  momentum  in  the  space  step  at  time  (t  )  to 

n+1 

known  values  in  the  space  step  from  the  previous  time  step  (t  ). 

r. 

Conservation  of  Energy 

To  obtain  the  finite  difference  form  of  the  Conservation  of  Energy 
equation,  begin  by  integrating  eq.  124. 


ftn+1  fXrm-1  [d  .  ,  U2  .  3  .  .  U: 

J  J  Jl  [p(e+2  )]  +  3^  ^u(e+2 


t  "  X 
n  m 


+*)  ]  dx  dt  =  0 
P 


rXm+1  rCn+1  d_  u2 .  i  ftn+1  1  3  u2  £  , 

I  J  3t[p(e+2  )]  dt  ax  +  I  J  3^[pu(e+2  +p)j  dx  dt 

JxJt  JtJx 

m  n  n  m 


Again,  eqs.  204  and  208  are  used,  with  the  following  result. 


,  ,  u2  .  m+1/2 

[/p (e+-  )] 


r  /  u  \ 

[ P ( e+2  ) , 


T  r  u2  p 

K[fRU(E+I  +i^i 


U2  P 

.RU(E+-  +i)]m 


Equation  211  relates  the  energy  in  the  space  step  at  time  (t  ) to  known 

n+ 1 

values  in  the  space  step  from  the  previous  time  step  (t  ). 

n 

Equations  209,  210,  and  211  form  the  basis  of  the  finite  difference 
technique  utilized.  Recall  from  figure  6  that  the  lower  case  letters 
indicate  space  averages  at  a  particular  instant,  and  the  upper  case  lctte: 
indicate  time  averages  at  a  particular  location. 


The  same  results  can  be  obtained  more  formally  (but  less  physically)  by 
use  of  the  two-dimensional  Green's  Theorem,  which  relates  a  surface  integral 
to  a  contour  integral. 


17 


Chapter  3 


I 

COMPUTATION  OF  THE  RIEMANN  PROBLEM 


The  essence  of  Godunov's  method  lies  in  the  technique  used  to  estimate 
the  assumedly  constant  properties  on  the  cell  boundaries,  as  depicted  by 
upper  case  letters  (e.g.  F  in  figure  6).  Each  boundary  is  considered  as  a 
Riemann  Problem  (or  shock  tube  analogy)  as  discussed  in  the  Introduction 
above.  The  equations  describing  the  constant  (in  time)  fluid  properties 
behind  the  waves  resulting  from  the  resolution  of  the  discontinuities  at  the 
cell  boundaries  will  now  be  derived.  The  two  types  of  waves  that  may  result 
(expansion  and  compression,  or  shock)  will  be  considered  separately. 

Shock  Wave  Equations 

Recall  that  each  computational  cell  is  assigned  an  average  value  for 
each  property,  and  this  average  value  is  considered  constant  throughout  the 
cell.  The  difference  in  value  from  one  cell  to  the  next  results  in  a  dis¬ 
continuity  at  the  common  cell  boundary.  If  the  resolution  of  the  discon¬ 
tinuity  (the  rupturing  of  the  diaphragm  in  the  shock  tube  analogy)  results 
in  a  shock  wave,  its  properties  may  be  computed  using  the  integral  equations 
of  motion,  equations  125,  126,  and  127,  as  shown  below,  (ref  5  and  ref  6). 

Note  that  the  three  integral  equations  of  motion  involve  integrals  of 
the  following  form. 

pa  j(  t ) 

I(t)  -  f (x , t )  dx 

Ja0(t) 

(301) 


where  f(x,t)  represents  p,  pu,  or  p(e+u2/2)  ,  and  is  discontinuous  at  x  =  , 

as  shown  in  figure  7.  Note  also  that  the  subscript  1  denotes  conditions 
ahead  of  the  discontinuity,  while  subscript  0  denotes  those  behind.  The 
conservation  laws  require  evaluation  of  the  derivative  of  eq.  301. 


r  <- ' 

\  •-j 


!(t> 

f(x,t) 


dx 


(302) 


Since  f(x,t)  is  discontinuous  at  x  =*■  the  integral  is  split  as  follows. 


£ 

dt 


|(t)-  € 

f (x, t) 

a0(t) 


dx  +  f  (x,  t)  [  (f+e )  -  (£ -O  ] 


I 


a1(t) 

f (x , t)  dx 
f  (  t  )  +€ 
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where  f(x,c)  is  (by  the  mean  value  theorem)  some  value  between  the  minimum 
value  of  f(x,t)  (fx)  and  the  maximum  value  (f0)  in  the  interval  [£±e],  That 
is,  f(x,t)  is  finite.  Since  the  term  in  the  large  curly  brackets  above  is 
a  limit  of  the  sum  of  three  well  behaved  functions,  the  order  of  the 
operations  may  be  reversed  and  the  sum  of  three  limits  may  be  substituted. 


-  tt 


■£<tW 


f ( x , t )  dx 


an(t) 


+  w 


2  €  f (x , t ) 


+  w 


(•a  L(t) 

f ( x , t )  dx 
£<t)+£ 


Since  f(x,t)  is  finite,  the  second  limit  term  is  zero.  Furthermore, 


r  r£ct)-e  -]  r<?C t> 

lira  f(x,t)  dx  -  f(x,t)  dx 

L  Jan(t)  J 


pa  ^  t ) 

j  f ( x , t )  dx 


'a0(t> 
r*a  !<  t ) 


r  J  f(x- 

J  J£(t) 


t)  dx 


— I(t)  -  — 
dt  ^  ;  dt 


written: 

•£(t) 

pa  t ) 

f(x, t) 

dx  +  dE 

^(t> 

a0ct) 

f(x,t)  dx 


(303) 


Leibnitz'  rule  (ref  3)  states  that  if 


pa  t ) 

J(t)  -  f ( x , t )  dx 

Ja0(t) 


then 


d_ 

dt 


pa  t(t) 

f(x,t) 

Ja0(t) 


dx 


pa  t ) 
Ja0(t) 


df(x, t) 

at 


dx  +  f(at 


f(a0,t) 


o 

dt 
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With  this,  eq.  303  becomes 


r^CO 
^  a0(t) 

pa^t) 

•'l(t) 


3f,(x,  t) 

at 


dx  +  f0(£,t) 


d£ 

dt 


af(x.t) 

at 


dx  +  f ( a  j , 


0 


where  f0(£(t),t)  and  f1(C<t)>t),  or  more  simply  f0  and  f  x  ,  are  shown 
in  figure  7.  Also, 


da0 

dt 


"  u0 


dat 

dt 


ui 


d£ 

dt 


W 


Note  that  W  is  the  wave  velocity,  measured  relative  to  the  fixed  frame  of 
reference  (the  x  axis),  and  make  the  following  substitutions. 


r€<t) 

Ja0Ct) 


af(x,t) 

at 


pa^t) 

dx  + 

J£<t) 


af(x.t) 

at 


dx 


+  f0W  -  f(a0(t),t)u  +  f(a1(t),t)u1  -  f ! W 


Now  shrink  the  interval  so  that  it  contains  only  the  discontinuity.  That 


let  a0 

-*  |  and  at  -» 

f .  Note  the  following  limits 

lim  | 

r?(t)  af(x.t) 

,  at 

dx  -  0 

lim 

ar«  J 

rai(t)  df(x.t) 
<««> st 

dx  -  0 

lim 

acT? 

f  (a0(t) ,  t)  j 

“  fo 

lim 

a  r« 

f  (at(t) ,  t)  ] 

J 

-  fi 

With  these  substitutions,  the  integral  equation  is  reduced  to  an  algebraic 
equation . 


^I(x,t)  -  0  +  0  +  f0W  -  f0u0  +  f i u j  -  fxW 
^I(x,t)  -  f  x  (ux  -W)  -  f o  (u0  -W) 


Combine  this  with  eq,  302  to  obtain  the  following  result. 


d_ 

dt 


r 


Lct) 


f (x, t) 


a0(t) 


dx  -  f x (ux -W) 


f0(u0-W) 


( 304) 


Recall  that  subscript  1  denotes  conditions  on  the  leading  edge  of  the 
shock  wave,  and  0  denotes  those  on  the  trailing  edge.  W  is  the  velocity  of 
the  shock  wave  relative  to  the  fixed  frame.  Equation  304  may  be  used  to 
reduce  the  integral  equations  of  motion  (eq's.  125,  126,  and  127)  to 
algebraic  equations  relating  the  flow  parameters  across  the  discontinuity. 

First  combine  eq's.  125  and  304  to  obtain  the  continuity  equation  in 
algebraic  form. 


Pi(ux-W)  -  p0(u0-W)  -  0 


Define 

w  -  W  -  ux 

(305) 

as  the  velocity  of  the  wave  relative  to  the  fluid  into  which  the  wave  is 
propagating.  With  this  substitution,  the  continuity  equation  becomes 

Pi[u1-(w+u1)]  -  p0 [u0-(w+ux) ]  -0 


■  1 

1 

“ 

+ 

ux  -  u0 

.  Pi 

P  0  . 

Piw  [  ~p  ]  +  [  u  ]  “  0 


(306) 


The  square  brackets  in  eq.  306  indicate  the  change 
quantity  across  the  shock  wave .  For  the  case  of  a 
(W  is  positive) ,  let  the  mass  velocity  of  the  wave 
over  by  the  wave  in  a  unit  time  per  unit  area  [ref 
where 


in  a  discontinuous 
positive  running  wave 
(the  mass  of  fluid  swept 
7])  be  denoted  by  b, 


b  “  P  y  W 

so  that  eq.  306  becomes 


(307) 


t [;]*[»]■« 


(308) 
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For  the  case  of  a  negative  running  wave  (W  is  negative) ,  the  mass  velocity 
is  denoted  by  a,  as  follows. 

a  -  -Pi w 

(309) 

Substituting  into  eq.  306  yields  the  conservation  of  mass  for  the  negative 
running  wave . 


a  [7]  '  [u]  -0 


(310) 


Combining  eq.  126  and  eq.  304  yields  the  momentum  equation  in  algebraic 

form. 


p iUi(ul-W)  -  p0u0(u0-W)  -  Po'Pi 


Substitute  eq.  305. 


PlUl [UX- (W+Ui) ]  -  PoUotUo-Cw+Ui)]  -  (P0-P1) 


Rearrange : 


-Pi wCUi-Uo)  +  (p  x "Po ) 


' 

Pi  W 

'  1 

1 

+  (uru0) 

L. 

.  Pi 

Po  . 

J 

The  quantity  in  curly  brackets  is  zero  by  the  continuity  equation  (eq.  306). 
P iw(Uj  *u0 )  -  (p 1 -p0 )  -  0 

(311) 

For  a  positive  running  wave,  substitute  eq.  307. 

b(ux -u0)  -  (Pi'Po)  “  0 
or  with  the  bracket  notation, 
b[  u  ]  -  [  p  ]  -  0 

(312) 

For  a  negative  running  wave,  substitute  eq.  309. 

-a(ux -u0)  -  (pi'Po)  -  0 


or  with  the  bracket  notation, 
a[  u  ]  +  [  p  ]  -  0 
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(313) 


Finally,  combining  eq.  127  and  eq.  304  yields  the  energy  equation  in 
algebraic  form. 


P i(ei+21)(ui'w)  -Po(eo+2°)(uo‘w)  ”  (Pouo-Piui) 


For  simplicity,  let 


e  + 


and  substitute  eq.  305. 


Pi'i  [u^Cw+Ui)]  -  p0£otuo-(w+ui)]  -  (Pouo'Piui) 


Rearrange : 


plw(e1-e0) 


(Piui-Pouo) 


r 

“  P o£o  • 

/>lW 

'  1 

1 

[ 

.  P\ 

Po  . 

(uj  -u0) 


} 


Again  the  term  in  the  curly  brackets  is  zero  by  the  continuity  equation, 
eq.  306. 


Piw(£i'£o)  -  (Px^i-PoUq)  -  0 

(314) 

For  a  positive  running  wave,  substitute  eq.  307.  With  the  square  bracket 
notation,  the  energy  equation  is  written 

b[e  +  J"]  -  [  pu  ]  -  0 

(315) 


For  a  negative  running  wave,  substitute  eq.  307.  With  the  square  bracket 
notation,  the  energy  equation  is  written 


a 


[  e  +  2“  ]  +  [  PU  ] 


-  0 


(316) 


The  shock  wave  equations  are  summarized  in  Table  1.  The  square  brackets 
denote  jump  conditions  across  the  wave,  and  now  subscripts  are  as  denoted  in 
figure  2. 
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Positive  Running  Wave 
b[i]  +  tu]  -0 
b  [  u  ]  -  [  p  ]  -  0 
b[e  +  2"]  -  [  pu  ]  -  0 
b  -  p4 w 


Negative  Running  Wave 

«[;]  -  [»]  -o 

a[u]  +  [p]“° 

aCe  +  2~^  +  Cpu]  “  0 

a  -  -prw 


Table  1.  Equations  relating  flow  parameters  across  a  discontinuity. 


The  next  task  is  to  evaluate  a  and  b.  To  do  so,  the  positive  running 
wave  and  the  negative  running  wave  must  be  considered  separately.  A 
negative  running  compression  wave  (figure  2)  is  characterized  by  p0  >  p4 . 
Equations  310  and  313  are  combined  to  eliminate  the  velocity  by  solving 
Eq.  310  for  [u]  and  substituting  the  result  into  eq.  313. 


a[a(J)]  +  [p]  -  0 


a 


P4-P0 


1/2 


(317) 


For  an  ideal  gas  (see  Appendix  A), 

1  _  (7-l)Po  _ +jLx±H&4  1  I 
P 3  L  (T+Dpo  +  (7-DP4  J  P\ 


When  this  is  substituted  into  eq.  317  and  the  terms  rearranged,  the  result 
is  the  following  equation  for  the  mass  velocity  of  the  negative  running 
compression  wave. 


a  “  V7P4P4 


(318) 


The  positive  running  compression  wave  (figure  2)  is  characterized  by 
Po  >Pi •  The  above  steps  are  repeated,  this  time  using  equations  308  and 
312.  Equation  308  is  solved  for  [u]  and  the  result  is  substituted  into 
eq.  312,  When  the  actual  pressures  and  densities  from  figure  2  are  used  in 
place  of  the  bracket  notation,  the  result  is  the  following. 
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Pi-Po  V2 

'1  -  1 
~P  1  ~P  2. 


(319) 


For  an  ideal  gas  (see  Appendix  A) , 

i  _  r  laillB I  +  (.7+1)  Po  1  1 
Pi  L  (7+l)Pi  +  (7-l)Po  J  P 2 


When  this  is  substituted  into  eq .  319  and  the  terms  rearranged,  the  result 
is  the  following  equation  for  the  mass  velocity  of  the  positive  running 
compression  wave. 


b  “  /yPi P\ 


(320) 


Expansion  Wave  Equations 

The  differential  equations  of  motion  are  valid  across  an  expansion  wave. 
To  reiterate,  these  are: 


k  <<■> + 1  (<,u>  - 0 


(113) 


3u  +  u  in  _  _  1 

at  ax  p  ax 


(119) 


/>(e  +  ?>] +  k  h(e  +  f +  ?>]  - 0 


(124) 


Equation  124  will  not  be  used  in  this  form.  It  will  be  replaced  by  the 
isentropic  gas  law, 


n 

7  —  constant 

P 


since  flow  across  an  expansion  wave  is  adiabatic  and  reversible. 
321,  p-p(p)  only,  so 


(321) 


From  eq . 


dv  _  dp  dp 
dx  dp  3x 


(322) 


Also,  one  can  show  (e.g.  ref  4)  that  for  weak  waves,  u  =  u (p)  only,  so 
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do  _  d o  du 
dt  ~  du  dt 

and 

dp  _  dp  du 
dx  ”  du  dx 

so  that  eq.  322  becomes 

dx>  dp  do  du 
dx  dp  du  dx 

Substitute  eq.  325  into  eq .  119. 

du  du  _  1  dp  dp  du 
dt  U  dx  p  dp  du  dx 


(323) 

(324) 

(325) 

(326) 


Now  expand  eq.  113, 


¥  *  -  0 


dt  dx  dx 


do  1  dp  dp  du  du 

du  p  dp  du  dx  ^dx 


But  0,  so  the  quantity  in  curly  brackets  must  be  zero.  That  term  is 

set  equal  to  zero  and  rearranged,  with  the  following  result. 
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Where  the  known  conditions  ahead  of  the  expansion  wave  (subscript  1  in  the 
case  of  a  positive  running  wave  as  shown  in  figure  2)  have  been  used  to 
evaluate  the  constant.  Solving  for  p  and  substituting  into  eq .  328, 


Now  integrate  over  the  whole  expansion  wave  in  order  to  relate  the 
conditions  behind  it  (subscript  0)  to  those  ahead  (subscript  1) . 
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Rearrange  this  and  substitute  the  isentropic  gas  law  (eq.  321)  in  the 
following  form: 


£■ o  _ 

P 1 


to  obtain 


Ui-U0 


+ 


2 

7-1 


’  £i" 

(t  1  )/2 1 

.  Poj 

When  the  right  hand  side  is  multiplied  by 


Ls.iz2.ol 

(Pi-Po) 


and  the  terms  rearranged,  the  equation  reduces  to  the  following  simplified 
form, 


±  K(Uj  -u0)  -  (p^Pq)  -  0 


(329) 


where 


K 


J 7P i P i 


■  1- 

P°/P1 

1- 

Kp,1 

L 

|<7-1>/27  ^ 

J 

(330) 


Again,  the  positive  sign  in  eq.  329  is  associated  with  the  positive  running 
wave  and  the  negative  sign  is  associated  with  the  negative  running  wave. 

From  figure  2,  conditions  ahead  of  the  positive  running  wave  are  denoted  by 
subscript  1,  while  conditions  ahead  of  the  negative  running  wave  are  denoted 
by  replacing  the  subscript  1  in  eqs .  329  and  330  with  a  subscript  4.  The 
pressure  and  velocity  behind  either  wave  is  denoted  by  subscript  0.  Using 
the  square  bracket  notation  to  represent  conditions  across  the  expansion 
wave,  eq.  329  can  be  rewritten 


±  K[u]  -  [p]  -  0 

(331) 


Finally,  let  a  -  -K  and  b  -  +K.  This  is  done  so  that  the  equations  are 
similar  to  the  momentum  equations  for  the  shock  wave,  eqs.  312  and  313. 
This  substitution  will  simplify  the  final  computer  program. 
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For  the  negative  running  expansion  wave  (p4>p0), 


a[u]  +  [p]  “  0 


where 


1- 

a  “  VtpTpT  ~ 


Po 


/P  4 


P°/P4j<7'1)/27 


and  for  the  positive  running  expansion  wave  (Pi>p0). 
b[u]  -  [p]  -  0 


where 


TtpIpT 


Po 


/Pi 


1- 


P°/Pi 


(7-D/27 


(332) 


(333) 


(334) 


(335) 


Generalized  Conditions  Behind  a  Wave 

Attention  is  now  turned  to  the  practical  problem  of  using  the  equations 
derived  above  to  solve  the  Riemann  problem  at  a  particular  cell  boundary. 
Notice  that  the  momentum  equation  for  the  negative  running  expansion  wave, 
eq.  332,  is  identical  in  form  to  that  of  the  negative  running  compression 
wave,  eq.  313.  The  only  difference  is  in  the  equation  used  to  evaluate  the 
mass  velocity  a.  Thus  any  numerical  scheme  used  will  have  to  test  to  see  if 
p4  is  greater  or  less  than  p0  and  then  use  the  appropriate  equation  to 
compute  a.  A  similar  arguement  can  be  made  for  the  positive  running  waves 
and  the  value  of  b.  With  these  differences  in  mind,  the  different  forms  of 
the  variables  a  and  b  are  momentarily  disregarded.  Thus  for  the  negative 
running  wave,  expansion  or  compression, 

a[u]  +  [p]  -  0 


or 

a(u4-u0)  +  (p4-p0)  -  0 


(336) 


and  for  the  positive  running  wave,  expansion  or  compression, 
b[u]  -  [p]  -  0 


or 

b (uj -u0 )  -  (Pi-Po)  ”  0 


(  337) 
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Solve  eq.  336  for  p0 . 

Po  -  au4  +  p4  -  au0 

(338) 

Solve  eq.  337  for  u0 . 

u0  “  k  (Po-Pi)  +  ui 

D  (339) 


Now  substitute  eq.  339  into  eq.  338. 

Po  -  au4  +  P4  -  §  (Po-Pi>  *  aui 

After  rearranging, 

bp4  +  apj  +  ab(u4-u1) 


Now  substitute  eq .  340  into  eq.  339  and  rearrange  to  obtain  the  following. 


uo 


au4  +  bu!  +  (p4-pi) 
a  +  b 


(341) 


Equations  340  and  341,  when  combined  with  the  equations  for  a  and  b, 
are  highly  nonlinear.  They  must  be  solved  iteratively  for  p0 ,  a  and  b, 
after  which  u0  is  found  easily.  The  iterative  technique  used  in  this  work 
is  discussed  in  detail  in  Appendix  B. 

In  the  case  of  weak  (sonic)  waves,  approximate  formulae  may  be  used  to 
compute  p0  and  u0 .  These  approximate  formulae  eliminate  the  need  for 
iteration  and  significantly  reduce  computer  time.  Experience  has  shown  that 
the  full  nonlinear  equations  need  only  be  used  when  the  cell  is  in  proximity 
to  an  actual  shock  wave  in  the  flow  field.  A  comparison  of  the  magnitudes 
of  pressure  in  the  two  adjacent  cells  suffices  to  determine  the  need  for  the 
nonlinear  equations  when  computing  the  Riemann  problem  on  any  particular 
boundary. 

Weak  waves  result  from  the  resolution  of  the  discontinuity  at  a  cell 
boundary  when  the  fluid  conditions  (velocity,  density,  pressure)  in  the  left 
and  right  cells  are  nearly  equal.  In  most  cases  the  grid  cell  size  can  be 
chosen  small  enough  that  the  average  properties  in  any  two  adjacent  cells 
are  nearly  equal.  In  this  way,  the  approximate  formulae  may  be  used 
everywhere  except  in  the  neighborhood  of  a  discontinuity  in  the  actual 
flowfield.  In  this  case  large  differences  in  properties  will  exist  across 
boundaries  of  adjacent  cells.  Note  that  the  distinction  between  expansion 
and  compression  waves  is  unimportant  for  weak  waves. 

The  essence  of  the  approximation  lies  in  the  assumption  that  weak  waves 
will  travel  at  the  local  speed  of  sound.  Thus  the  velocity  of  the  negative 

running  wave  is 
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w  -  Jrpjfu 


In  terms  of  the  mass  velocity  a  (eq.  309)  this  is 


w  -  a/p4 


Equating  the  two  forms  yields 


a  -  y-rp 4P< 


(342) 


for  a  sonic  wave.  The  pressure  and  density  in  the  negative-most  cell  can  be 
written 


P4  +  Pi  P4  *  Pi 
P4  “  - 0 -  +  - 0 - 


P  4  +  P 1  P  4  ‘  P 1 


For  small  differences  in  pressure  and  density,  the  second  term  in  each 
expression  is  seen  to  be  much  smaller  than  the  first,  and  may  be  neglected. 
Substituting  the  remaining  expressions  into  eq.  342  gives  the  mass  ve loci  tv 
of  the  negative  running  sonic  wave. 


f  P4  +  Pi  1 

•  7  [  2  . 


P  4  +  Pi 


(343) 


For  the  positive  running  wave,  the  mass  velocity  is 

b  -  Jy p ~ p j 

where  the  pressure  and  density  in  the  positive-most  cell  may  be  written  as 

Pi  +  P4  Pi  '  P4 
Pi  -  O  +  - 0 - 


Pi  +  P 4  Pi  '  P 4 


Again  Che  second  terms  are  neglected,  and  the  remainders  used  in  the 
expression  for  b. 


f  P4  +  Pi  -]  P*  +  Pi  -]  V' 

b  ”  '  7  [  2  J  L  2  J  ’ 


(344) 


So  for  small  differences  in  the  fluid  properties  in  the  two  cells,  the  mass 
velocities  of  the  two  waves  are  considered  equal.  This  is  valid  because  the 
density  and  pressure  are  nearly  equal  in  the  two  cells.  Combining  eq.  343 
and  eq.  344, 


a  =  b 


{"N 1  - 

■  P4  +  Pi  [  p*  +  Pi  1 

L  — J| 


(345) 


With  a  -  b,  the  momentum  equations  for  the  negative  and  positive 
running  weak  waves  are : 


aH  +  [p]  ”  0 


(negative  running  wave) 


aH  -  [p]  -  0 


(positive  running  wave) 


With  the  subscript  notation  of  figure  2,  these  are  written 
a(u4-u0)  +  (p4-p0)  -  0 


(346) 


a(u1-u0)  -  (Pi-Po)  -  0 


(347) 


Equation  347  is  subtracted  from  eq.  346  and  the  result  rearranged  to  obtain 


Pi  +  P< 


Equation  347  is  added  to  eq .  346  and  the  result  rearranged  to  obtair 


(348) 


u,  +  u4 


P 4  -  Pi 


+  a  2 


(349) 


In  summary,  the  linear  approximations  to  the  solution  of  the  Riemann  problem 
at  a  cell  boundary  are: 


Pi  +  P  4 

Po  “  - 2  +  a 


2 
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(348) 


ui  +  U4  iP«  ■  Pi 
2  +  a  2 


a 


7 


P4  +  Pi  1 
2 


Pi  +  P in 
2 


(349) 


(345) 


These  equations  apply  to  weak  waves,  whether  they  are  expansion  or 
compression,  negative  running  or  positive  running.  Unfortunately,  the 
situation  is  slightly  more  complicated  for  determination  of  density.  As 
seen  in  figure  2,  the  density  behind  the  waves  is  different  on  each  side  of 
the  contact  discontinuity.  The  density  behind  the  negative  running  wave 
(expansion  or  compression)  is  (see  Appendix  A) 


P  3 


(7+l)Pn  +  (.7-1)  Pi 
.  (7*l)Po  +  (7+l)P4  . 


P  4 


(350) 


and  the  density  behind  the  positive  running  wave  (expansion  or  compression) 
is  (see  Appendix  A) 


X7±1.).Eo..+..L7-U.Bi 
.  (7-l)Po  +  (7+l)Pi 


Pi 


(351) 


The  Computational  Scheme 

A  typical  Riemann  situation  is  shown  in  figure  8.  The  discontinuity  at  cell 
boundary  m  has  resolved  itself  into  two  waves.  The  conditions  ahead  of  the 
waves  are  known  from  the  previous  time  step.  The  conditions  at  the  cell 
boundary  (denoted  by  uppercase  letters)  are  sought,  so  that  fluxes  across 
this  boundary  may  be  determined.  Four  possible  outcomes  exist,  depending  on 
the  direction  of  propagation  of  each  wave.  The  fluid  properties  behind  the 
wave  are  given  in  Table  2  for  each  case.  In  all  cases,  the  energy  is  found 
from  the  equation  of  state.  Using  the  subscript  notation  of  figure  2,  and  a 
and  b  for  the  mass  velocities,  the  absolute  velocities  of  the  waves  are 


w 


u, 


Pi 


(positive  running  wave) 


(352) 


w  -  u4 


a 

Pi 


(negative  running  wave) 


(353) 
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Case  1.  w+  and  w  are  both  positive. 


U  -  u4 

R  -  ,4 

P  “  P4 


Case  4.  w+  and  w  have  different  signs;  u0 


U  -  u0 
R  -  P2 
P  -  Po 

is  negative . 


Table  2.  Fluid  properties  behind  the  waves  emanating  from  the  cell 
boundary  for  all  four  possible  conditions. 
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Chapter  4 


SAMPLE  PROBLEM  1 

Shock  Wave  Reflection  from  a  Wall 


Theoretical  Study 

The  reflection  problem  may  be  stated  as  follows.  Fluid  flows  steadily 
in  a  one  dimensional  pipe  at  a  velocity  of  Uj  *  V,  as  shown  in  figure  9a. 

At  x  —  0  and  t  —  0 ,  an  ideal  valve  is  instantaneously  closed,  stopping  the 
fluid  motion  at  that  point.  The  problem  is  to  determine  the  resulting  fluid 
motion  (figure  9b) . 

Nature  reconciles  this  seemingly  impossible  situation  by  creating  a 
compressive  wave  in  the  incoming  fluid.  This  wave  travels  into  the  incoming 
fluid,  and  is  precisely  of  the  proper  strength  required  to  leave  the  fluid 
behind  itself  motionless  (the  same  velocity  as  the  wall). 

In  deriving  the  equations  describing  this  process,  a  positive  sense  is 
assumed  for  both  the  incoming  fluid  and  the  reflected  wave,  even  though  they 
must  be  of  opposite  sense  in  any  real  problem.  As  with  all  derivations,  the 
proper  technique  is  to  assume  positive  for  the  derivation,  and  substitute 
positive  or  negative  quantities  (as  the  case  may  require)  into  the  resulting 
equations.  The  benefit  of  this  in  the  present  case  is  to  allow  a  single 
derivation  for  both  reflection  problems:  the  wall  at  the  left  (incoming 
fluid  velocity  is  negative  and  reflected  wave  velocity  is  positive)  and  the 
wall  at  the  right  (incoming  fluid  velocity  is  positive  and  reflected  wave 
velocity  is  negative) . 

The  following  notation  applies  to  the  flow  field.  The  variables  are 
shown  schematically  in  figure  9b. 

Subscript  1:  conditions  ahead  of  the  wave 
Subscript  0:  conditions  behind  the  wave 

u:  absolute  velocity  of  the  fluid 

W:  absolute  velocity  of  the  shock  wave 

v:  velocity  of  the  fluid  relative  to  the  wave 

Equation  305,  relating  fluid  velocity  relative  to  the  wave,  will  be  required 
in  the  following  forms . 


v0  -  u0  -  W 


(400) 


vi  -  u,  -  W 


(401) 
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From  Appendix  A,  the  conservation  equations  relating  conditions  across  a 
shock  wave  are : 


PlVj.  -  PoV0 

Pi  +  Pivi  -  Po  +  P ovo 

e  +  a  o  "f"  +  ^Vq 

Pi  Po 

Also  from  Appendix  A,  the  equation  of  state  of  a  perfect  gas  is: 


(501) 

(502) 

(503) 


e 


_1_  E 
7-1  p 


From  figure  9b,  ux  -  V  and  u0  »  0 ,  so  the 
the  wave  (eqs.  400  and  401)  are  as  follows. 


fluid  velocities 


(511) 

relative  to 


v0  -  -W 

(402) 

vx  -  V  -  W 

(403) 

The  conservation  equations  (eqs.  501  through  503)  can  be  reduced  to  a 
more  useable  form,  as  follows.  Begin  by  substituting  eqs.  402  and  403  into 
eq.  501  and  rearranging. 


Po 


Pi  (1 


(404) 


Substitute  eq.  402  and  eq .  403  into  eq .  502  and  rearrange,  then  substitute 
eq.  404  to  obtain  the  following. 


Po  “  Pi  +  Piv2(l  -  f) 


(405) 


To  modify  the  energy  equation  (eq.  503),  substitute  eqs.  402  and  403.  The 
equation  of  state  (eq.  511)  is  written  twice  (once  for  subscript  0  and  once 
for  subscript  1)  and  substituted  into  eq.  503.  Next,  eqs.  404  and  405  are 
substituted,  with  the  following  result. 


2~L  Ei 
7*1  L  Pi 


+  ^..WV  +V2  -2WV  -  0 
7-1 


Finally,  divide  through  by  V,  since  V  *  0  (trivial  solution),  multiply  by 
(V-W) ,  and  rearrange  to  obtain  the  following  quadratic  equation  in  W. 
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The  solution  of  the  quadratic  equation  is 


W  -  V  ~  [  V)2  +  V2  +  j 

(406) 

For  the  present,  wave  reflection  from  a  wail  at  the  left  is  to  be  computed. 
For  this  case,  the  incoming  velocity  is  negative,  and  we  are  looking  for  a 
positive  wave  velocity.  For  this  example  problem,  the  following  values  are 
given  (ref  8) : 

V  -  -428  m/s  (to  the  left) 

Pl  -  1.225  kg/m3 
Pl  -  101,300  n/m2 
7  ~  1.4 

When  these  values  are  substituted  into  eq.  406,  the  wave  velocity  is  either 
of  the  following. 


W 


f  +255  m/s 
\  -590  m/s 


Consistent  with  the  problem  statement,  the  positive  velocity  is  selected. 
This  is  substituted  into  eqs .  404  and  405  to  obtain  the  remaining  properties 
of  the  reflected  wave. 

W  -  255  m/s 

p0  -  3.280  kg/ra3 

p0  -  459,400  n/m2 


Numerical  Study 

The  reflection  problem  stated  above  was  solved  by  the  Godunov  method  in 
order  to  test  its  ability  to  predict  shock  wave  reflection  from  a  wall. 

This  phenomenon  plays  an  important  part  in  the  problem  which  ultimately  will 
be  solved;  a  hypervelocity  jet  travelling  through  a  cylinder. 

The  one  dimensional  computational  grid  is  shown  in  figure  10.  It  shows 
the  left  and  right  boundaries,  and  a  typical  cell  boundary  (m) .  Each  fluid 
property  in  the  domain  (velocity,  density,  pressure,  energy)  is  discretized 
by  assigning  some  average  value  in  each  cell.  These  values  are  constant 
throughout  the  cell,  and  are  denoted  by  lower  case  letters  (u,p,p,e).  As  a 
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result  of  this  discretization,  a  discontinuity  exists  at  each  cell  boundary. 
The  discontinuity  is  resolved  by  expansion  or  compression  waves  that  emanate 
from  the  boundary.  The  constant  fluid  properties  behind  each  set  of  waves 
(i.e.  at  the  boundary)  are  denoted  by  upper  case  letters  (U,R,P,E).  Sub¬ 
scripts  for  average  cell  properties  (u,p,p,e)  are  denoted  by  the  letter  i, 
while  subscripts  for  boundary  values  (U,R,P,E)  are  denoted  by  the  letter  m. 

Boundary  conditions  are  provided  by  virtual  cells  (shown  dashed  in 
figure  10)  at  the  wall  and  upstream.  Reflection  conditions  at  the  wall  (m  - 
0)  are  provided  by  assigning  density,  pressure,  and  energy  in  the  virtual 
cell  to  be  equal  to  those  of  the  first  cell,  while  the  velocity  in  the 
virtual  cell  is  set  to  the  negative  of  the  velocity  in  the  first  cell.  This 
ant isymme trie  condition  guarantees  that  the  fluid  velocity  at  the  wall  will 
be  zero. 

On  the  downstream  boundary,  free  stream  conditions  are  specified  in  the 
virtual  cell  for  all  time.  This  boundary  condition  may  cause  a  disturbance 
to  propagate  into  the  computational  domain.  Consequently,  the  domain  was 
made  large  enough  to  insure  that  the  disturbances  would  not  reach  the 
computational  region  of  interest  for  the  length  of  time  being  computed. 

The  following  initial  values  are  provided  to  begin  the  computation. 

The  energy  was  found  from  the  equation  of  state,  eq.  511. 

u.  -  -428  m/s 

p  -  1.225  kg/m3 

p.  -  101,300  n/m2 

The  flow  chart  for  the  one  dimensional  code  is  shown  in  figure  11.  A 
listing  of  the  code  is  given  in  Appendix  C  (this  is  actually  the  code  for 
sample  problem  2,  the  shock  tube,  so  the  boundary  conditions  are  different). 

Results  are  presented  graphically  in  figures  12  through  15.  The  code 
predictions  are  in  excellent  agreement  with  the  theoretical  determinations 
made  above.  Figure  12  shows  the  fluid  velocity  at  three  different  times 
after  reflection.  At  t  -  .001  s  for  example,  the  reflected  wave  is 
about  0.25  m  from  the  wall.  The  fluid  velocity  ahead  of  the  wave  is  428  m/s 
(to  the  left),  while  the  fluid  velocity  behind  the  wave  is  zero.  The  wave 
propagates  to  the  right  at  a  velocity  of  256  m/s  (the  theoretical  value  was 
255  m/s).  The  shock  is  spread  over  one  cell  (0.10  m) ;  the  fluid  velocity  is 
zero  at  x  -  0.20  m  and  -428  m/s  at  x  -  0.30  m.  This  typical  result  is  the 
best  possible  resolution  with  this  finite  grid. 

The  density  shock  wave  is  seen  in  figure  13.  The  predicted  density  in 
the  first  cell  is  low  by  about  5%,  but  all  other  cells  are  within  1%.  The 
problem  in  the  first  cell  is  discussed  by  Godunov  in  ref  2  (paragraph  6) . 

He  believes  that  the  system  "errs  in  entropy"  in  computing  the  unsteady 
portion  of  the  reflection  process.  No  attempt  was  made  here  to  correct  the 
problem,  since  it  only  affects  the  first  cell.  The  pressure  shock  (figure 
14)  and  the  energy  shock  (figure  15)  are  also  in  excellent  agreement  with 
theory  (within  1%),  although  the  energy  curve  shows  a  slight  distortion  at 
the  wall  due  to  the  density  error  there. 
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Chapter  5 


SAMPLE  PROBLEM  2 
The  Shock  Tube 


Experimental  Results 

An  experiment  was  conducted  in  the  BRL  2 -inch  (51  mm)  shock  tube  to 
check  the  predictions  of  the  one -dimensional  code.  Figure  16  is  a 
photograph  of  the  shock  tube  and  the  instrumentation  used.  The  main 
components  of  the  instrumentation  shown  in  figure  16  are  quartz  pressure 
transducers  and  a  four  channel  digitizing  oscilloscope.  Dimensions  of  the 
shock  tube  are  given  in  figure  17.  The  ends  of  the  shock  tube  were  closed 
in  order  to  reflect  the  waves. 

Each  test  began  with  air  at  atmospheric  pressure  and  temperature  in  each 
section  of  the  shock  tube.  The  two  sections  were  separated  by  a  0.002  inch 
(0.05  mm)  mylar  diaphragm.  Pressurized  air  at  room  temperature  was 
introduced  slowly  into  the  driver  secton  until  the  diaphragm  ruptured. 

Figure  18  shows  a  typical  diaphragm  after  rupture.  A  total  of  five  tests 
were  done,  but  the  results  of  the  first  were  erratic  so  were  disregarded. 

The  data  from  the  remaining  four  tests  was  averaged  and  this  average  was 
used  to  compare  to  the  code  prediction. 

The  table  below  gives  the  data  that  was  manually  recorded  for  each  test, 
and  the  averge  value  used  in  the  code  verification. 


Test 

2 

3 

4 

5 

AVG 

Driver  Pressure,  psig 

61.5 

59.1 

60.1 

58.6 

59.8 

Room  Temperature,  °C 

25.3 

26.4 

27.1 

27.7 

26.6 

Atmospheric  Pressure,  psia 

14.83 

14.83 

14.83 

14.83 

14.83 

At23,  ns 

212. 

215. 

213. 

216. 

214. 

The  output  from  the  pressure  transducers  was  stored  in  the  digitizing 
oscilloscope.  Immediately  following  each  test,  the  data  was  reduced  by  a 
microcomputer  (shown  in  figure  19)  and  stored  on  tape  for  later  plotting. 
These  plots  are  given  in  figures  20  -  23,  the  pressure  -  time  trace  for  each 
gage  on  each  test.  These  data  were  also  averaged,  and  plotted  as  the 
circles  in  the  graphs  at  the  end  of  this  section. 

The  utilization  of  the  data  in  figures  20  -  23  begins  in  time  with  the 
arrival  of  the  first  shock  wave.  This  is  of  necessity  t  -  0,  since  the 
oscilloscope  does  not  trigger  when  the  diaphragm  breaks,  but  rather  writes 
continuously.  Similarly,  data  after  about  8  msec  (station  1)  or  9  msec 
(station  2)  is  not  useful,  as  it  is  altered  by  a  small  shock  wave  that  is 
reflected  from  the  remains  of  the  diaphragm  (resulting  from  the  main 
reflected  wave  which  is  by  this  time  travelling  back  into  the  driver  section 
after  its  reflection  from  the  wall  in  the  driven  section) .  Data  from 
station  3  was  not  used  in  this  evaluation. 
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Experimental  Chamber  Pressure  Correction 

While  the  code  assumes  an  ideal  diaphragm  rupture,  in  reality  it  never 
is.  Some  of  the  energy  of  the  gas  in  the  driver  (chamber)  is  consumed  in 
this  non- ideal  rupture.  This  pressure  loss  is  well  known  to  shock  tube 
researchers,  having  been  discussed  in  1950  by  C.  W.  Lampson  (ref  9)  of  BRL. 
This  loss  can  be  corrected  when  predicting  shock  pressure  by  reducing  the 
measured  chamber  pressure  by  some  amount,  and  then  assuming  a  perfect 
diaphragm  rupture.  To  determine  this  "diaphragm  opening  coefficient", 
theoretical  shock  pressures  must  be  computed.  (Other  factors  contributing 
to  inaccuracies  are  viscous  effects,  including  the  effect  of  the  boundary 
layer  on  the  gage  reading,  and  the  finite  volume  of  gas  in  the  driver 
section. ) 


To  obtain  an  independent  comparison,  the  following  nonlinear  equation 
for  shock  pressure  is  taken  from  Anderson  (ref  4) . 


Pi 


1  - 


(7-l)(}1)<S°-l) 

a4  Pi 

27[27+(7+D(^-l)]] 


•,-27/7-1 


where  p4  -  initial  absolute  chamber  pressure 

p,  -  initial  absolute  driven  section  pressure 

p0  -  pressure  behind  the  resultant  shock  wave 

Since  the  gas  in  both  sections  of  the  shock  tube  is  initially  at  the  same 
temperature,  the  sonic  velocities  are  equal:  at  -  a4 .  For  air  (7  =  1.4)  the 
equation  reduces  to 

p  (Po/Pi)  -  1  y7 

P4  “  Po  1  ~  ~  ~  .. 

L  y7[7+6((Po/Pl)-l)]J 


For  the  problem  at  hand,  pt  is  atmospheric  pressure  (.102  MPa  during  this 
experiment) . 

To  determine  a  theoretical  curve  of  shock  pressure  (p0)  vs  chamber 
pressure  (p4),  one  needn't  solve  the  above  non-linear  equation  for  p0 . 
Instead,  values  of  p0  are  substituted  and  p4  is  computed.  The  results  of 
this  are  plotted  in  figure  24.  Also  shown  is  experimental  data  supplied  by 
G.  Coulter  of  TBD,  BRL.  The  difference  is  apparent.  Close  examination  of 
the  curves  shows  that  the  experimentally  determined  chamber  pressure  should 
be  reduced  by  5.5%  (in  the  region  .480  <  p  <  .540  mpa)  to  account  for  losses 
in  the  non-ideal  diaphragm  rupture.  That  is,  in  an  experimentally 
determined  chamber  pressure  of  p4 ,  5.5%  of  it  is  consumed  by  the  diaphragm 
rupture,  leaving  .945  p4  to  form  the  shock  wave.  Consequently  the  chamber 
pressure  used  in  the  code  must  be  .945  p4 ,  where  p4  is  the  experimental 
chamber  pressure,  and  .945  is  the  diaphragm  opening  coefficient. 
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Theoretical  Study 


Definition  of  Terns.  The  following  schematics  illustrate  the  subscript 
notation  used.  As  can  be  seen,  4  denotes  the  driver  (high  pressure)  side, 
and  1  the  driven  side. 


nit;ial 
ltions: 


Selection: 


ftSJt 


ef feet f on: 


U4 

P4 '  p 

>  e4 

U1 >Pl 

-ei 

diaphragm 

4 

- 

- 1 

!  u0,p0 

1 

2 

1 

1 

exp  wave 

contact  surface 

incident  shock 

- 1 - 

1 

1 

u0,Po 

2 

w  - 

5 

1 

—  » 

r 

reflected  shock 


Data.  The  following  measurements  were  taken  during  the  experiment. 

T4  -  26.6°C  Tj  -  26.6°C 

p4  -  74.6  psia  (.487  MPa)  p*  -  14.8  psia  (.102  MPa) 

In  addition,  the  initial  velocities  in  the  two  chambers  are  assumed  to  be 
zero . 

u4  -  0  Uj  -  0 


Computed  Conditions.  The  following  values  are  required  for  the  theoretical 
computations  and  as  inputs  to  the  computer  code.  A  value  of  287  J/kg*K  was 
used  for  R,  and  of  course  7  -  1.4. 


P 4  “  ”  5,66  k§/m3 

P\  ”  -  1.19  kg/m3 

e<  -  ^  RT4  -  215,100  J/kg 
e 1  -  ^71  RTi  “  215,100  J/kg 
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Theoretical  Computations 


In  order  to  avoid  any  biasing  of  the  theoretical  predictions  in  favor  of 
the  code  predictions,  Anderson's  equations  (ref  4)  will  be  used  instead  of 
the  relations  derived  in  this  report  for  the  one  dimensional  code. 

A.  Speed  of  Sound  in  the  Undisturbed  Gas ,  a1  and  a4 . 

at  -  VtRTY  -  347  m/s 
a4  -  -  347  m/s 

B.  Pressure  of  gas  behind  the  incident  shock  wave,  p0  (ref  4,  eq.  7.94). 


£4  _  Eo  J 
Pi  Pi 


1  "  T 


(7-l)(T1)(S°-l) 

a4  Pi 


27[27+(7+l)(?°-l)] 
Pi 


-27/tI 


Substitute  7  -  1.4  and  ax  «  a4  and  rearrange: 

0.4(p0-p1) 


Po  “  P4 


1  - 


y2 . 8p  x  72.4po-fO.4pj 


After  substituting  p,  and  p4 ,  the  above  equation  is  solved  by  iteration  with 
the  result: 

p0-  .213  MPa 

C.  Velocity  of  incident  shock,  W.  (ref  4,  eq  7.14). 

W.  -  [  (Eo-l)  +  1 

>  L  27  vpx  J 

W.  -  482  m/s 


D.  Velocity  of  gas  behind  the  incident  shock  wave,  u0  (ref  4,  eq.  7.16). 

i* 


u0  -  (*>-l) 

7  Pi 


_2 ^ 
_i±L 


Eo  +  Xii 
Pi  7+1 


u0  -  194  m/s 
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E.  Density  behind  the  incident  suock  wave  ( p2 )  (ref  4,  eq.  7.11). 


P  2  ”  P  1 


1  +  0±i  £0 

- 1-iPi- 

Xti  20 

T-l  Pi 


p2  -  1.991  kg/m3 


F.  Pressure  behind  the  reflected  shock  wave  (p5).  The  velocity  of  the 
incoming  fluid,  V  —  u0 ,  is  positive  in  the  current  problem.  Equation  406  is 
now  applied  to  a  reflected  wave  travelling  in  the  negative  x  direction  (to 
the  left),  thus  requiring  that  W  be  negative.  Equations  404  and  405  are 
then  used  to  compute  the  density  and  pressure  behind  the  reflected  wave.  In 
the  current  notation,  eq.  406  is  written  as  follows. 


W  -  -326.6  m/s  (to  the  left) 

W 

Ps  "  Po  +  P2v2(l  "  ^ ) 


ps  -  .414  MPa 


(405) 


G.  Density  behind  reflected  wave  (ps). 


r 


(404) 


P5  -  3.174  kg/m3 


Numerical  Study 

The  computational  grid  used  for  the  shock  tube  problem  is  presented 
schematically  in  figure  25.  A  complete  listing  of  the  code  used  to  solve 
this  problem  is  given  in  Appendix  C.  The  code  results  are  presented 
graphically  in  figures  26  and  27,  and  numerically  in  the  summary  of  the 
experimental  evaluation  that  is  given  in  Table  2  below. 

Figures  26  and  27  show  the  pressure  -  time  history  at  stations  1  and  2 
respectively.  The  solid  line  shows  the  code  prediction,  while  the  circles 
indicate  experimental  data.  The  code  prediction  is  based  on  the  corrected 
chamber  pressure  in  both  figures.  The  circles  represent  averages  of  the 
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data  from  the  four  tests  (figures  20  -  23) .  Since  the  instrumentation  did 
not  record  the  diaphragm  rupture  (i.e.  t  -  0) ,  the  experimental  arrival  time 
of  the  incident  shock  was  made  to  coincide  with  that  predicted  by  the  code. 

In  figure  26,  the  incident  shock  wave  is  seen  to  arrive  at  station  1  at 
1.2  msec  after  diaphragm  rupture.  The  constant  pressure  behind  the  shock 
wave  begins  to  drop  at  about  2.5  msec,  as  the  expansion  wave  (reflected  from 
the  driver  end  of  the  shock  tube)  arrives.  At  4.5  msec,  the  reflected  shock 
(reflected  from  the  closed  end  of  the  driven  section)  passes  through  the 
expansion  wave  at  station  1,  causing  a  jump  in  pressure.  The  expansion  wave 
continues  to  reduce  the  pressure  at  station  1,  when  at  5.2  msec  the  head  of 
the  expansion  wave  returns  (after  reflection  at  the  closed  end  of  the  driven 
section  of  the  shock  tube) ,  causing  a  noticeably  increased  rate  of  pressure 
decrease.  At  about  8.5  msec  the  experimental  traces  showed  a  reflection  of 
the  reflected  shock  wave  from  the  remains  of  the  diaphragm,  so  no  more 
experimental  results  are  plotted.  The  code,  however,  shows  the  main  driver 
end  rf -  lection  of  the  reflected  shock  wave,  which  arrives  at  9.3  msec. 


In  figure  27,  the  incident  shock  wave  is  seen  to  arrive  at  2.4  msec, 
followed  quickly  by  the  reflected  shock  (2.9  msec),  and  the  expansion  wave 
(that  was  reflected  from  the  driver  end)  at  3.6  msec. 

Both  figures  show  general  agreement  between  the  code  and  the  experiment. 
The  pressure  behind  the  incident  wave  is  predicted  by  the  code  to  be  a 
constant  ,213  MPa  at  both  stations,  while  it  actually  was  .215  MPa  at 
station  1  and  .212  MPa  at  station  2.  The  code  predicts  the  pressure  behind 
the  reflected  wave  to  be  .413  MPa,  while  it  actually  was  .407  MPa. 

Collision  of  waves  is  similarly  well  predicted,  as  can  be  seen  in  figures  26 
and  27.  Verification  of  this  fact  was  a  principal  goal  of  this  study. 


Digital  output  from  the  code  was  used  to  compute  the  velocities  of  the 

incident  (w  )  and  reflected  (w  )  shock  waves  between  stations  1  and  2. 

'  l  r 


w.  - 
1 

,1684m  - 

0.6096m 

00237s  - 

.00122s  ” 

w  - 

r 

.6906m  - 

1.1684m 

00446s  - 

,00290s 

486  m/s 
-358  m/s 


(to  the  left) 


The  incident  wave  speed  was  experimentally  determined  by  measuring  with  an 
electronic  counter  the  wave's  time  of  passage  between  stations  2  and  3. 


W.  - 


.2700m  -  1.1684m 
214/is 


-  475  m/s 


No  direct  measurements  were  made  for  the  reflected  wave  speed,  but  it  can  be 
obtained  from  the  pressure-time  traces  in  figures  20-23.  The  average 
velocity  of  the  reflected  wave  is 


_ 0.6906m  -  1.1684m _ 

.00122+. 00332)s  -  (. 00237+ . 00052 ) s 


-339  m/s 
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The  expansion  wave  properties  were  not  predicted  as  well  as  the  shock 
wave  properties  and  the  wave  velocities.  This  is  common  in  shock  tube 
predictions  (for  example,  see  ref  10  and  11)  and  is  believed  to  be  the 
result  of  the  interaction  of  the  cold  (expanded)  gas  from  the  driver  section 
with  the  hot  (compressed)  gas  in  the  driven  section.  Although  the  code  does 
not  consider  this  effect,  the  wave  shapes  in  figures  26  and  27  are  accurate, 
and  the  magnitudes  of  the  pressure  predictions  are  within  10%,  and  usually 
much  less  than  that. 

Note  that  a  small  error  in  arrival  time  of  the  reflected  wave  (at 
station  1  for  example)  causes  a  large  error  in  the  prediction  of  expansion 
wave  pressures.  In  fact  some  of  the  error  in  the  expansion  wave  predictions 
seen  in  figures  26  and  27  may  be  due  to  viscous  effects  in  the  boundary 
layer  of  the  expansion  wave.  Any  slowing  of  the  wave  would  cause  the 
experimental  points  in  the  expansion  wave  to  be  shifted  to  the  right  with 
respect  to  the  code  prediction. 

The  code  predicts  shock  wave  speeds  to  be  slightly  higher  than  the 
experimentally  determined  values  with  the  wave  travelling  the  farthest  (the 
reflected  wave)  being  most  in  error.  Also,  a  drop  in  the  experimental 
pressure  behind  the  incident  wave  is  noticed  between  stations  1  and  2 .  Both 
of  these  phenomena  are  attributed  to  viscous  attenuation  of  the  shock  wave . 
If  viscous  effects  were  included  in  the  code,  the  predictions  of  wave  speed 
and  pressure  would  be  even  better  than  those  discussed  above.  However, 
viscous  effects  are  rightfully  negligible  in  the  hypervelocity  problem 
ultimately  to  be  solved,  so  no  attempt  was  made  to  modify  the  code.  Some 
quantification  is  possible,  however,  and  is  useful  to  the  overall  evaluation 
of  the  code. 

Emrich  and  Wheeler  (ref  12)  discuss  shock  attenuation  due  to  wall 
effects  (mainly  a  turbulent  boundary  layer)  and  propose  the  following 
empirical  equation  to  predict  the  shock  attenuation. 


Eo  .  i 

wm 

1 - 

1 

o 

C4 

L _ 

.  Pi 

X 

L  Pi  J 

exp(-A^) 


where  D  -  hydraulic  diameter  of  the  shock  tube 
x  -  distance  the  shock  has  travelled 
A  -  an  empirical  constant 
p0  -  pressure  behind  the  shock  wave 
px  -  pressure  ahead  of  the  shock  wave 

They  give  examples  of  the  use  of  this  equation  for  shock  strengths  (p0/Pi) 
of  1  to  15.  They  also  suggest  a  value  of  0.0024  for  A.  The  shock  strength 
in  the  experiment  under  study  here  is  about  5.  To  evaluate  the  constant  A, 
measurements  of  p0  for  the  incident  wave  made  at  stations  1  and  2  are  used. 


station 

x/D 

Eo 

Ei 

1 

12 

.215  MPa 

.  102  MPa 

2 

23 

.212 

.  102 
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These  values  are  substituted  into  the  equation,  using  (23-12)  -  11  diameters 
for  the  distance  the  shock  has  travelled.  The  result  is  that  A  =  0.00245. 
Thus  the  general  relationship  is 


_Eo_  . 

.102 


_Ro_ 

.102 


exp(-. 00245  |  ) 


This  method  reconciles  the  difference  between  the  experimental  reflected 
pressure  at  station  2  (.407  MPa)  and  the  code  prediction  (.413  MPa),  which 
was  close  to  the  theoretical  prediction  (.414  MPa).  However,  the 
experimental  value  for  the  pressure  behind  the  incident  shock  at  station  1 
(.215  MPa)  was  actually  higher  than  theoretical  and  code  predictions  (both 
were  .213  MPa).  This  is  contrary  to  the  effects  of  viscous  attenuation  and 
must  be  attributed  to  experimental  error  (which  includes  the  possibi  lity  of 
inputting  erroneous  initial  conditions  into  the  code) . 

Attenuation  of  the  expansion  wave  and  its  reflections  was  not 
considered,  due  to  uncertainty  of  the  applicability  of  this  method  to  the 
pressure  distribution  within  the  expansion  wave. 

The  experimental  verification  of  the  one -dimensional  code  is  summarized 
in  table  2.  The  summary  shows  that  the  code  predictions  are  in  excellent 
agreement  with  both  theory  and  experiment  for  all  flow  parameters  studied. 

The  code  predictions  were  also  checked  against  the  predictions  of  a 
Beam-Warming  type  of  code  used  in  the  Terminal  Ballistics  Division  of  the 
BRL  to  study  and  predict  shock  tube  behavior.  The  BRL  uses  shock  tubes 
extensively  in  the  study  of  explosive  blast  dynamics;  how  shock  waves  from 
explosions  affect  structures,  for  example. 

For  this  comparison,  code  predictions  for  velocity,  pressure  and  density 
as  functions  of  position  in  the  shock  tube  were  plotted  every  millisecond 
from  0  to  12  milleseconds .  These  plots  are  shown  in  figues  23,  29  and  30, 
where  the  plot  for  t  -  12  milleseconds  is  omitted  due  to  space  restrictions . 
Looking  at  the  pressure  plots  (figure  30)  for  example,  one  can  see  the 
original  shock  wave  travelling  to  the  right,  and  the  original  expansion  wave 
travelling  to  the  left  at  t  —  1  millesecond.  The  expansion  wave  has  already 
begun  its  reflection  from  the  left  wall.  At  t  -  3  milleseconds  the  shock 
wave  has  reflected  from  the  right  wall.  The  two  reflected  waves  collide 
between  3  and  4  milleseconds,  followed  by  reflection  of  the  expansion  wave 
from  the  right  wall  beginning  just  prior  to  t  -  4  milleseconds.  The 
reflected  portion  of  the  expansion  wave  interacts  with  the  incident  portion 
in  the  t  -  .005  millesecond  frame,  and  begins  to  chase  the  reflected  shock 
wave  that  is  travelling  toward  the  left  wall.  The  two  actually  collide  just 
prior  to  impact  with  the  left  wall.  The  remaining  frames  show  a 
continuation  of  this  process  of  reflection  and  collision. 

One  additional  observation  worthy  of  note  is  seen  in  the  density  plots 
(figure  29).  The  contact  surface  is  clearly  seen  at  t  *  1  millesecond 
between  the  left  running  expansion  wave  and  the  right  running  shock  wave. 
This  contact  surface  moves  to  the  right,  and  its  interaction  with  the 
reflected  shock  wave  is  seen  in  the  t  -  4  millesecond  and  t  -  5  millesecond 
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Experimental 

Code 

Theoretical 

Pressure  behind  incident 

Shock  wave  (p0).  MPa 
(station  1/station  2) 

.215/. 212 

.213 

.213 

Pressure  behind  reflected 

Shock  wave  (ps)  MPa 

.407 

.413 

.414 

Speed  of  incident 

Shock  wave  (W.  )  m/s 

475 

486 

482 

Speed  of  reflected 

shock  wave  (W  )  m/s 
r 

339 

358 

327 

Velocity  of  gas  behind  the 
incident  wave  (u0)  m/s 

- 

193 

194 

Velocity  of  gas  behind  the 
reflected  wave  (u,.)  m/s 

- 

<0.01 

0 

Density  of  gas  behind  the 
incident  wave  p2  kg/m3 

- 

1.99 

1.99 

Density  of  gas  behind  the 
reflected  wave  ps  kg/m3 

- 

3.16 

3.17 

Table  3.  Summary  of  results  showing  a  comparison  between  code  predictions, 
theoretical  predictions,  and  experimental  data. 


frames.  Following  the  collision,  the  contact  surface  moves  to  the  left. 
Notice  that  the  code  predicts  a  smearing  out  (or  smearing)  of  this  contact 
surface  with  time. 

Figure  31  shows  the  comparison  of  the  predictions  of  the  two  codes. 
Velocity,  density  and  pressure  are  plotted  as  functions  of  position  for 
times  of  4,  8  and  12  milleseconds .  These  are  the  same  as  the  plots  just 
discussed.  Predictions  of  the  Godunov  type  code  developed  here  are  again 
shown  as  a  solid  line.  Superimposed  on  these  plots  are  the  results  of  the 
Beam-Warming  code,  shown  as  circles.  In  every  case  the  agreement  between 
the  two  codes  is  excellent. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


This  effort  was  the  first  phase  of  a  project  to  write  an  axisymmetric 
fluid  dynamics  code  and  apply  it  to  the  hypervelocity  problem  of  a  shaped 
charge  jet  within  a  cylindrical  tunnel.  In  this  phase  of  the  project,  the 
necessary  equations  were  derived  and  a  one  dimensional  Godunov  code  was 
written  for  an  inviscid,  perfect  gas. 

Theoretical  and  experimental  verifications  were  done  which  showed  that 
the  Godunov  technique  accurately  predicted  those  flow  phenomena  critical  to 
this  study:  shock  and  expansion  wave  formation,  wave  reflection,  and 
wave /wave  collisions.  The  details  are  summarized  as  follows: 

1.  Predictions  of  pressure  were  within  2%  of  the  experimental  values. 

2.  Predictions  of  velocities  were  within  6%  of  the  experimental  values 

3.  Predictions  of  fluid  densities  were  within  1%  of  theoretical  values 

As  a  result  of  the  techniques  learned  and  verifications  made,  a 
decision  was  made  to  expand  the  present  code  to  include  axisymmetric  and 
real  gas  considerations.  This  effort  is  reported  on  in  ref  21. 
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FIGURES 
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Figure  lb.  Discretized  Pressure  in  a  One  Dimensional  Fluid. 
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Figure  lc.  The  Classical  Shock  Tube  Problem. 


Figure  1.  The  Discretized  Flowfield 
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Figure  2 .  The  x  -  t  Diagram  for  the  Riemann  Problem 


3 

I 


a 

i 


iLUfi.LLi. 


4J 


I 


s 


U 

C 

o 

o 

ai 


rrifiTTi 

+* 

X 

w 

a 


53 


Figure  4.  The  Lagrangian  Parameters 
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Figure  5.  A  Pressure  Discontinui 


Figure 
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Figure  7.  A  Snapshot 
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Figure  8.  A  Typical  Computation  Situation 


Figure  9a.  Before  reflection. 


Figure  9b.  After  reflection. 


Figure  9,  The  One  Dimensional  Reflection  Problem 
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SHOCK  WAVE  REFLECTION  FROM  A  WALL 


61 


Figure  12.  The  Velocity  Shock  Wave 


SHOCK  WAVE  REFLECTION  FROM  A  WALL 
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Figure  14  .  The  Pressure  Shock  Wave 


SHOCK  WAVE 
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Figure  15.  The  Energy  Shock  Wave 


Figure  16.  A  Photograph  of  the  Shock  Tube 
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Figure  17.  A  Schematic  of 


Figure  18  .  A  Ruptured  Diaphragm  from  the  Shock  Tube 
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Figure  24,  Pressure  Behind  the  Shock  Wave  in  a  Shock  Tube 
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Computational  Domain  for  the  Shock  Tube 
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i  26.  Pressure  at  Station  1 


Figure  27.  Pressure  at  Station  2 
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Figure  28.  Predicted  Fluid  Velocity  in  the  Shock  Tube 
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igure  29.  Predicted  Fluid  Density  in  the  Shock  Tube 
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Figure  30.  Predicted  Fluid  Pressure,  in  the  Shock  Tube 
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lines:  the  present  Code 

circles:  Beam-Warming  Code 


Position,  meters 
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Appendix  A 


DENSITY  RATIO  ACROSS  THE  WAVE 


In  deriving  the  density  ratio  across  a  wave,  the  shock  wave  will  be 
considered  first.  Following  that  derivation,  some  comments  will  be  made 
regarding  the  applicability  of  using  the  same  expression  for  computing  the 
density  ratio  across  an  expansion  wave. 

The  shock  wave  equations  for  a  negative  or  a  positive  running  wave  were 
derived  in  Chapter  3  (eqs.  306,  311,  and  314).  Those  equations  relate  the 
conditions  across  the  shock  wave,  and  do  so  in  terms  of  absolute  velocities 
of  the  fluid.  The  equations  may  be  expressed  in  terms  of  fluid  velocities 
relative  to  a  moving  shock  wave,  since  the  conditions  are  invariant  with 
respect  to  shock  wave  translation  (ref  13) .  This  transformation  is  detailed 
in  ref  4,  and  the  result  is  given  here.  The  following  notation  is  required. 

Subscript  1:  conditions  ahead  of  the  wave 
Subscript  0:  conditions  behind  the  wave 

u:  absolute  velocity  of  the  fluid 

W:  absolute  velocity  of  the  shock  wave 

w:  velocity  of  the  wave  relative  to  the  fluid  ahead 

When  expressed  in  terms  of  the  above  notation,  the  shock  wave  equations 
become  the  following. 


P iw  -  p 0(W-u0) 

Pi  +  p\y>2  =  p0  +  p0  (w-u0 ) 2 

ei  +  +  ^w2  =  e0  +  ^°  +  ^(W-u0)2 

Pi  P  o 


(501) 

(502) 

(503) 


Equations  501  through  503  relate  conditions  behind  a  moving  shock  wave 
(negative  or  positive  running)  to  conditions  ahead  of  the  shock  wave,  and  do 
so  in  terms  of  fluid  velocities  relative  to  the  moving  wave.  Notice  that 
(W-u0)  is  the  velocity  of  the  shock  wave  relative  to  the  fluid  behind,  or 
the  negative  of  the  velocity  of  the  fluid  behind  the  shock  wave  relative  to 
the  shock.  Similarly,  w  is  the  negative  of  the  velocity  of  the  fluid  ahead 
of  the  shock  wave  relative  to  the  shock. 


These  equations  are  manipulated  to  obtain  the  density  ratio  across  the 
shock  as  a  function  of  the  pressure  across  the  shock.  Equations  501  and  502 
are  used  to  obtain  expressions  for  the  velocities  w  and  (W-u0),  each  in 
terms  of  the  pressures  and  densities  of  the  fluid  ahead  of  and  behind  the 
shock  wave.  These  expressions  are  then  used  to  eliminate  the  velocities 
from  the  energy  equation,  eo  503,  with  the  following  result. 
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e,  +  El  +  4 

Rir — Ro 

=  e0  +  Eo  +  * 

Rir — Ro 

1  Pi  *  Pi 

.Pi"  Po_ 

P  0  P  0 

.Pi"  Po_ 

(508) 


The  last  step  is  to  express  the  fluid  energy  in  terms  of  pressure  and 
density,  leaving  only  pressures  and  densities  in  the  equation.  For  a 
calorically  perfect  gas, 


p  -  pRT 

e  =  c  T 
v 

where  T  is  the  absolute  temperature  in  Kelvins.  Equations  509  and  510 
together  with  the  definitions 

R  —  c  -  c 
P  v 

7  ”  c  /c 
P  v 

can  be  combined  to  obtain  the  equation  of  state  of  a  perfect  gas. 


(509) 

(510) 


e  -  -i-  E 

7-1  P  (51 

The  equation  of  state  is  the  closure  equation  for  the  system  of  equations 
eq.  501  through  eq.  503  (three  equations,  four  unknowns).  Equation  511  is 
substituted  into  eq.  508,  and  after  considerable  rearranging,  the  desired 
result  is  obtained. 


.  .  (7+l.lPo  ±  ( 7 *  1 ) P i 

P°  (7-l)Po  +  (7+DPi  Pl 


(512) 


Equation  512  relates  the  density  behind  the  shock  wave  to  the  density  ahead 
of  the  shock  wave  (presumeably  a  known  quantity) ,  in  terms  of  the  strength 
of  the  shock  wave  (indicated  by  the  pressures  on  either  side  of  the  shock 
wave ,  pj  and  p0) . 

Even  though  the  shock  wave  equations  were  used  in  deriving  eq .  512,  no 
assumptions  were  necessary  regarding  the  strength  of  the  shock.  Therefore, 
eq.  512  is  appropriate  for  any  strength  shock  wave. 

Equation  512  is  known  as  the  shock  adiabat  (ref  1,14).  The  correspond¬ 
ing  relation  for  an  expansion  wave  is  the  Poisson  adiabat  (ref  1,14),  or 
ise^ntropic  gas  law,  eq.  321.  It  can  be  written  in  a  form  similar  to 
eq.  512: 


P  o 


Pi 


(513) 
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In  the  accoustic  approximation  (i.e.  for  waves  in  which  p0  =  pj),  the  two 
expressions  (eq.  512  and  eq.  513)  are  approximately  equal,  since  the  density 
ratio  tends  to  unity  in  each  case  (as  should  be  expected).  Further,  since 
in  the  limit  of  an  infinitely  weak  wave  there  is  no  distinction  between 
expansion  and  compression  waves,  use  of  eq.  512  for  both  types  of  weak  waves 
is  appropriate.  Godunov  (ref  1)  makes  the  assumption  that  eq .  512  can  be 
used  also  for  strong  expansion  waves,  making  it  the  only  equation  required 
to  compute  the  density  ratio  across  any  wave.  Experience  here  indicates 
that  this  is  a  valid  approximation  for  a  perfect  gas.  Also,  Peyret  and 
Taylor  (ref  15)  list  an  alternate  expression  for  the  density  and  state  that 
its  results  are  consistent  with  eq.  512. 
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Appendix  B 

ITERATIVE  SOLUTION  TO  THE  RIEMANN  PROBLEM  EQUATIONS 


In  most  cases  a  basic  iterative  approach  to  the  solution  of  the  non¬ 
linear  equation  for  the  pressure  behind  the  waves  emanating  from  the  cell 
boundaries  is  adequate.  These  equations  are  derived  above,  and  are  repeated 
here  for  convenience. 


Po 


bp4  +  apj  +  ab(u4-u1) 
a  +  b 


(340) 


The  mass  velocities  of  the  two  waves  that  emanate  from  the  boundary  under 
study  are  a  and  b. 


a 


Xzl 
2  7 


JlPiPi 


1- 


(v  1  )/2-» 


PO  <  P4 


a 


JwTpI 


Po  ^  P< 


(333) 


(318) 


b 


JlPxPi 


1- 


(7-D/27 


Po  <  Pi 


(335) 


b  ”  ■J'lPxPi 


Po  *  Pi 


(320) 


For  such  nonlinear  algebraic  equations ,  of  the  form 


Po  “  f(Po) 


(601) 


convergence  of  the  iterations  is  not  always  guaranteed.  Godunov  (Ref  2) 
indicates  the  equation  may  not  converge  for  a  strong  rarefaction  wave,  but 
presents  no  supporting  data.  Considering  the  complexity  of  the  equation, 
his  statement  may  have  been  based  on  actual  numerical  examples.  In  view  of 
this,  an  investigation  of  the  conditions  for  convergence  of  eq.  601  is 


35 


enlightening,  and  will  lead  to  a  modification  to  eq.  601  which  will 
guarantee  convergence  in  all  cases,  and  in  general  will  speed  up  the  rate  of 
convergence.  A  summary  of  the  method  used  is  given  below,  followed  by  a 
detailed  derivation  of  the  method. 


Summary  of  the  Method 

The  condition  for  convergence  of  the  iterations  of  the  nonlinear 
algebraic  equation,  eq.  60 1,  is  expressed  mathematically  by 


f'(o!  <  i 


(610) 


where  f  is  the  solution  of  eq.  601.  Computing  the  derivative  of  eq .  340  is 
a  horrendous  task.  Fortunately,  modifications  can  be  made  to  eq.  601  which 
will  guarantee  convergence  of  any  function,  for  any  finite  value  of  f ' ,  and 
without  any  knowledge  of  the  analytic  form  of  f ' .  In  one  approach,  Chorin 
(ref  17)  and  Sod  (ref  18)  use  the  following  iterating  function  in  place  of 
(eq .  601) . 


Po  -  [l-c]p0  +  of(p0) 


(614) 


Notice  chat  the  p0  on  the  right  hand  side  represents  the  current  value  of 
the  iterate,  while  the  p0  on  the  left  hand  side  represents  the  new  value. 
Chorin  and  Sod  assign  arbitrary  values  to  the  constant  c  (but  c  <  1)  until 
convergence  is  achieved.  However,  convergence  of  eq.  614  is  guaranteed  if 
c  satisfies  the  following  equation. 


C  =  l-f'(f) 


(617) 


Substituting  eq.  617  into  eq .  614  results  in  the  following  iterating 
function. 


£iBo)-fitnp0 
i  -  f'(n 


(618) 


With  superscripts  to  denote  sequential  iterates,  the  derivative  is 
approximated  in  the  following  manner. 


f(Po  )  ~  fCPo1) 

i  i-1 

Po  Po 


1 ' 
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This  is  consistent  with  the  accuracy  of  the  finite  difference  method  that  is 
being  used  to  numerically  solve  the  fluid  equations  of  motion.  Now  eq .  618 
is  rewritten  with  iterative  superscripts  and  eq.  621  is  substituted,  and  the 
following  iteration  function  results. 


i+1 

Po  ‘ 


f(p’o  Mp'o  -  Po  Af5 
Ap'o  -  Af* 


(623) 


where 


i-l 

Po 


(624) 


Af*  -  f  (p'0  )  -  f  (p'o  1 ) 


(625) 


Equation  623  is  the 
computations,  p°  and  Po 


recommended  iterating  function.  To  begin 
are  required.  These  may  be  computed  with 


the 


Po 


Pi  +  P  4 
2 


(626) 


Po  -  f (Po) 


(627; 


Details  of  the  Derivation 

Derivation  of  the  conditions  for  convergence  of  eq .  601  can  be  found  in 
most  texts  on  numerical  analysis  (e.g.  ref  16),  ar.d  is  summarized  in  the 
following  way. 

The  solution  of  eq.  601  is  denoted  by  f .  Then  in  eq .  601, 


r  -  f(r) 


(602) 


In  wuat  will  be  referred  to  below  as  the  standard  iteration  technique,  the 
iterations  are  begun  with  an  initial  guess,  p§  .  The  first  revised  value 
is  the  result  of  substituting  p°  into  eq.  601,  and  is  written 


Po  *  f(Po) 


87 


The  revised  value  is  then  substituted  into  eq.  601  for  a  new  revised  value, 
written  and  so  on.  With  the  help  of  the  mean  value  theorem,  one  can  show 
that  the  error  in  the  nth  iteration  is  a  function  of  the  derivative  of 
eq.  601,  f ' ,  as  follows. 


[f-pn0]  -  tr-po  1  n  f'(fn ) 

0  n 

(609) 

Here  f  represents  some  value  between  the  true  root  (C)  and  the  initial 

guess  (p°).  Thus  if  the  derivative  of  f  in  eq.  601  is  less  than  one  in  the 

neighborhood  of  the  root  f,  the  right  hand  side  of  eq.  609  will  include  the 

product  of  a  series  of  n  numbers,  each  of  which  is  less  than  one.  Thus  the 

error  in  the  nth  iteration, 

[C-Pn0  1 


will  be  less  than  the  error  in  the  initial  guess, 


ir-pg] 


and  this  error  can  be  made  as  small  as  desired  by  increasing  n.  Notice  that 
the  initial  guess  must  be  within  the  region  near  the  root  f  in  which  the 
derivative  is  less  than  one.  Also  notice  that  the  smaller  (i.e.  closer  to 
zero)  that  f'  is,  the  faster  will  be  the  convergence. 

The  condition  for  convergence  of  the  iterations  of  the  nonlinear 
algebraic  equation,  eq.  601,  is  expressed  mathematically  by 


|f'(f)|  <  1 


(610) 


Computing  the  derivative  of  eq.  340  is  a  horrendous  task.  As  mentioned 
above,  Godunov  states  that  the  iterations  of  eq.  340  may  not  converge  for 
strong  rarefaction  waves.  Five  numerical  examples  of  strong  rarefaction 
waves  were  computed  as  part  of  this  study,  and  in  all  cases  standard 
iteration  of  eq.  340  resulted  in  convergence  in  fifteen  or  fewer  iterations. 
Other  investigators  (ref  14,17,18)  indicated  they  had  encountered 
convergence  problems.  Modifications  can  be  made  to  eq .  601  which  will 
guarantee  convergence  of  any  function,  for  any  finite  value  of  f ' ,  and 
without  any  knowledge  of  the  analytic  form  of  f ' .  Godunov  (ref  2)  proposes 
such  a  modification,  but  without  any  supporting  rationale.  The  origin  of 
his  scheme  is  obscure .  Although  it  apparently  works,  several  authors  (ref 
14,17,18)  have  used  alternate  methods.  The  method  presented  here  is  similar 
to  Godunov's,  and  is  an  extension  of  the  method  used  by  Chorin  (ref  17)  and 
Sod  (ref  18).  It  is  apparently  similar  to  one  proposed  by  Bedijn,  as 
referred  to  in  appendix  A  of  ref  14,  but  no  paper  by  Bedijn  on  the  subject 
could  be  found. 
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The  first  step  in  deriving  the  modified  iterating  function  is  to  rewrite 
eq .  601. 


^(Po>  "  Po  -  f(Po) 

(611) 

The  solution  to  eq.  611  is  p0  -  f,  for  which  f  -  f(f)  (eq.  602  above),  and 


iKD  -  0 

<  61  ?) 

Solve  eq.  611  for  f(p0)  and  take  the  derivative: 
f'  -  „  i  . 

dpo 


Substitute  this  into  eq.  610 


1  -  V>'(p0)|  <  1 


o  <  r  <2 


(613) 


This  is  the  condition  for  convergence  of  the  iterating  function,  eq.  611. 

It  is  the  same  as  the  condition  (eq,  610)  for  the  iterating  function  eq. 
601,  in  that  the  function  f  in  eq.  611  is  still  subject  to  eq.  610  in  order 
that  eq.  613  will  be  satisfied.  However,  0  in  eq.  611  can  be  multiplied  by 
some  non-zero  constant  c  without  changing  the  root  f  in  eq.  612.  If  c  is 
properly  chosen,  eq.  613  can  be  satisfied  regardless  of  the  value  of  f' ,  A 
new  iterating  function  is  written. 


F(Po)  “  Po  -  cV>(p0) 

-  Po  -  c[Po'f(P0) ] 


Notice  that  at  the  root  f ,  the  quantity  in  brackets  is  zero  and  F(f) 
converges  to  the  root  f .  Rearranging, 


F(Po)  "  [l-c]po  +  cf(p0) 


(614) 


This  is  the  iterating  function  used  by  Chorin  (ref  17)  and  Sod  (ref  18). 
More  will  be  said  about  their  application  later. 
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To  study  the  conditions  for  convergence  of  the  iterating  function  F,  eq. 
614  is  differentiated  and  substituted  into  eq.  610. 


F’  -  [1-c]  +  cf' 
f'  -  \  [F'-l+c] 


and  in  eq.  610: 

ij  [ F ' * 1+c ] )  <  1 


I  F '  !  <  1 

(615) 

Not  surprisingly,  the  requirement  for  convergence  is  the  same  as  that 
of  f(p0).  The  improvement  results  from  the  arbitrary  constant  c,  whose 
value  can  be  chosen  such  that  eq .  615  is  satisfied  regardless  of  the  value 
of  f ' .  An  arbitrary  requirement  is  now  imposed  to  define  the  arbitrary 
constant  c,  the t  in  the  neighborhood  of  the  root  f  , 


F’(n  -  0 


(616) 


Surely  this  will  satisfy  the  convergence  criterion  eq.  615.  It  will  in  fact 
provide  the  quickest  convergence  of  the  iterations.  differentiate  eq.  614, 


F' (Po )  -  [1-c]  +  cf'(p0) 


F'(C)  -  [1-c]  +  cf'(r)  -  o 


1 

c  -  l-f'(f) 


(617) 


Selection  of  c  in  accordance  with  eq.  617  will  guarantee  convergence  of 
F  in  eq.  614.  By  inspection  of  eq.  617,  two  restrictions  on  f'($")  are 
apparent,  but  both  are  easy  to  accept  without  limiting  the  procedure. 

First,  since  c  *  0 ,  f'(f)  must  be  finite.  Second,  since  c  must  be  finite, 
f'(C)  *  1.  There  is  negligible  chance  of  this  occurring  numerically  (i.e. 
f'  -  1.000000  etc.  or  even  close  to  it)  in  a  function  as  complex  as  eq.  340. 
A  special  algorithm  for  computing  c  can  be  invoked  if  f'  is  computed  to  be 
too  close  to  1  for  reliable  computation  of  c.  This  was  not  a  problem  in  the 
present  work.  In  any  case  the  value  of  c  can  be  controlled  (limited  to  some 
maximu"’  value,  for  examp1 e)  since  there  is  some  latitude  in  the  choice  of 
the  magnitude  of  F' :  it  need  only  be  less  than  one,  not  necessarily  equal 
to  zero. 
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The  following  cases  are  possible  for  eq.  617. 


i.  f'(r)  <  o 


0  <  c  <  1 


Underrelaxation:  eq.  614  produces  a  new  value  between  p’0  and  f(p’0). 


2.  f'(f)  -  0  -  -  c  -  1 

Equation  614  reduces  to  the  original  function,  which  converges  unaided. 

3.  0  <  f'  (f)  <  1  -  -  c  >  1 

Overrelaxation:  new  value  from  eq.  614  is  outside  interval  p’0  to  f(p'0) 

4.  f '({■)-  1  -*  -»  c  -*■  “ 

Problem  area:  see  above  discussion. 

5.  f '($■)>  1  -  -  c  <  0 

convergence  :  Slgn  of  (p°  '  po  >  1S  °PPoslte  that  01  (Po  *  Po  )  • 


Equation  617  can  be  substituted  into  eq.  614  to  obtain  a  single 
expression  for  the  iterating  function. 


Po  “  F(p0) 


1  -  f'(f) 


(618) 


Note  that  combining  eq.  611  with  eq.  618  leads  to  the  Newton-Raphson  method 
(see,  for  example,  ref  19).  For  clarity,  let 


a  -  f'(f) 
so  that 


(619) 


Po  -  F(p0) 


HEq)  -ap0 
1  -  a 


(620) 


This  is  the  iterating  function  to  be  used  in  place  of  eq .  601.  To  verify 
that  the  solution  of  eq .  601  is  also  a  solution  of  eq .  620,  substitute 
Po  **  f(p0)  “  f  into  eq.  620. 


-  r 
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For  a  final  check,  differentiate  eq.  620,  where  a  -  £'(£)  is  not  a  function 
of  the  variable  p0 . 


F'(Po) 


f.llBo)  -  -  -Q 
1  -  a 


and  at  the  root  p0  —  f : 


F'(0 


f '(C)  -  a  a  -  a 
1  -  a  1  -  a 


Thus  the  iterating  function  F(p0)  will  converge  rapidly. 

Godunov  (ref  2)  proposes  an  iterating  function  similar  to  eq.  620,  but 
with  a  replaced  by  -a.  The  equation  he  gives  for  a  bears  little  analytical 
resemblance  to  -f'(S’),  and  r-°  numerical  resemblance.  The  derivation  of  his 
scheme  must  differ  from  the  above  derivation  in  some  way  (perhaps  many 
ways),  but  since  he  makes  no  comments  about  the  origin  of  his  scheme,  and 
since  the  origin  is  by  no  means  obvious,  it  must  remain  unknown  for  the 
present.  Thus  a  comparison  of  the  method  presented  here  with  Godunov's 
method  is  not  possible.  A  discussion  of  the  present  technique  for  improving 
convergence  of  the  iteration  may  be  found  in  ref  20. 

Chorin  (ref  17)  proposed  a  method  based  on  eq .  614.  Sod  (ref  18)  also 
used  the  method.  In  Chorin' s  method,  eq .  614  is  the  iterating  function,  and 
c  is  initially  assigned  a  value  of  1.  Since  this  reduces  eq.  614  to  the 
original  iterating  function  (eq.  601),  Chorin  starts  all  iterations  with  the 
original  function.  If  this  function  fails  to  converge  after  twenty 
iterations,  he  resets  c  to  1/2  and  continues  to  iterate.  If  convergence 
still  does  not  occur  in  another  twenty  steps,  he  resets  c  to  1/4,  and  so  on. 
In  practice,  he  states  he  never  had  to  go  beyond  1/4.  Thus  he  is  manually 
imposing  values  for  c  that  correspond  to  f'(f)  <  0  (cases  1  and  2  above,  for 
eq.  617).  Theoretical  justification  of  this  approach  is  not  discussed  in  his 
paper,  and  is  not  easily  obtained  (if  at  all)  from  340.  Numerically, 
however,  the  method  apparently  worked  well  for  him  as  well  as  for  Sod. 

The  difficulty  with  the  method  proposed  above,  and  that  which  is  avoided 
in  Chorin's  method,  is  that  the  value  of  a  =  £'  (f)  is  unknown,  due  to  the 
complexity  of  f(p0)  (eq.  340)  and  to  the  fact  that  the  solution  (f)  is  not 
known  a  priori.  But  since  the  iterations  are  simply  a  series  of  approx¬ 
imations  to  the  root  f,  using  a  valid  approximation  to  the  derivative  £'(?) 
is  quite  consistent.  Thus  eq.  619  is  written 


a 


f'(D 


df 

dPo 


Po”f 


Af 

Ap0 


Po=r 


i  _Af  1 
Ap'o 


f(Po  ) 


-  fCPo'1) 

i  - 1 
Po 


(621) 
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Now  rewrite  eq.  620  with  iterative  superscripts: 


?+1  f(Po  )  -  ®  (Po) 

1  -  a 

Substitute  eq.  621  into  eq.  622  to  obtain 
i+1  f(p'o  Mp'o  -  Po  Af5 

P0  A  ' 

Ap0  -  Af 

where 

»  >  1  i  •  1 

Ap  o  ”  P  o  "  V o 

Af3'  -  f(Po  )  -  fCpV1) 


(622) 


(623) 


(624) 

(625.'' 


Equation  623  is  the  recommended  iterating  function.  To  begin  the 
computations,  pg  and  pj  are  required.  These  may  be  computed  with 


Po 


Pi  +  P< 


(626) 


Po  -  f(Po) 


(627) 


Equation  626  is  suggested  by  several  authors  as  a  good  initial  estimate.  px 
and  p4  are  the  pressures  in  the  right  and  left  cells  prior  to  resolution  of 
the  discontinuity  at  their  common  boundary.  Equation  627  is  simply  the 
first  revised  estimate  based  on  the  standard  (unmodified)  iterating 
technique  (eq.  601).  Thus  the  first  two  iterates  are  the  same  as  those  of 
both  the  standard  method  (eq.  601)  and  Chorin's  method  (eq.  614). 

Van  Leer  (ref  14)  notes  that  during  the  iteration  large  negative  values 
of  [u4-ut]  (i.e.  ux  »  u4)  may  lead  to  a  negative  value  for  p0 .  Hence  he 

suggests  a  minimum  allowable  (positive)  value  for  p’0  ]  say  e,  which  is  the 
minimum  meaningful  pressure  for  the  flowfield  under  study.  If  the  next 
iteration  also  results  in  a  negative  pressure,  he  suggests  stopping  the 
iteration  and  setting  p0  -  e . 


Chorin  makes  another  numerical  suggestion.  He  notes  that  erroneous 
convergence  occurs  immediately  if  pj  -  p4  in  eq.  626.  Therefore  he  suggests 
that  his  iteration  scheme  be  carried  out  for  at  least  two  steps.  This 
requirement  is  guaranteed  in  the  method  presented  here,  since  the  first 
iteration  computed  by  equation  eq .  623  is  pj; .  This  is  the  first  iterate  to 


93 


be  compared  to  its  previous  iterate  to  determine  if  the  desired  level  of 
convergence  has  been  achieved. 

Several  numerical  examples  were  studied  using  the  four  methods 
discussed  above:  Unmodified  iteration,  Godunov's  method,  Chorin’s  method, 
and  the  present  method.  The  results  in  all  cases  substantiated  the 
conclusions  drawn  here  regarding  the  actual  convergence  of  the  various 
methods  for  different  situations,  as  well  as  the  rates  of  convergence.  In 
all  cases,  the  present  method  converged,  and  did  so  the  fastest.  The  other 
methods  did  not  converge  in  every  case. 

Two  of  the  numerical  examples  are  worth  noting.  The  first  was  simply  an 
algebraic  equation  in  which  f'  was  clearly  greater  than  one.  This  example 
demonstrated  non-convergence  of  both  the  standard  technique  (eq.  601)  and 
Chorin's  method  (eq.  614),  and  rapid  convergence  for  the  method  presented 
here  (eq.  623).  The  second  numerical  example  studied  a  strong  rarefaction, 
and  all  four  methods  (including  Godunov's)  converged.  The  method  presented 
here  converged  fastest. 
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Appendix  C 

LISTING  OF  THE  COMPUTER  PROGRAM 


PROGRAM  DIAFRAM 
C 

C  *********************************************************************** 
C 

C  10  APR  1985 

C  EXPERIMENTAL  VERIFICATION  OF  ONE-D  CODE 

C  HUBERT  W.  MEYER  JR. 

C 

c  ********************************************************************** 

c 

LOGICAL  L 

INTEGER  T , TT , STA1 , STA2 

PARAMETER  (MD-102 , MR-423 ,MT-525 , TT-2800 , STA1-305 , STA2-491) 
DIMENSION  CU(0:MT,0: 1) , CR ( 0 : MT , 0 : 1 ) , CP(0 :MT, 0 : 1) , CE(0 :MT , 0 : 1) 
DIMENSION  U(0:MT+I,0-1) ,R(0 :MT+I , 0 : 1) , P(0:MT+1 , 0 : 1) , E(0 :MT+1 , 0 : 1) 
DIMENSION  TIME(0:TT) 

DIMENSION  PCD  (0:50),FPCD  (0:50) 

REAL  MOM1 , MOM2 , MOM3 , NRG 1 , NRG 2 , NRG 3 
C 

T-0 

TIME(O)— 0 . 

H-.003 
PMIN-1 . 

G-1.4 

G1-(G-1)/(2.*G) 

G2-(G+1)/(2.*G) 

C 

C  GENERAL  NOTES: 

C  1.  IN  THE  DIMENSIONED  VARIABLES  U(I,T),  ETC,  THE  CURRENT 

C  TIME  STEP  IS  DENOTED  BY  T-0  AND  THE  NEW  (NEXT)  TIME  STEP  BY  T-l.  THUS 
C  U(17 ,0)  IS  THE  CURRENT  (KNOWN)  VALUE  AND  U(17,l)  IS  THE  NEW  VALUE 
C  BEING  COMPUTED.  THIS  NOTATION  APPLIES  TO  EACH  TIME  STEP  IN  TURN  AS 
C  THE  SOLUTION  IS  MARCHED  OUT. 

C  2.  IN  THE  PROGRAM,  I  ALWAYS  DENOTES  CELL  NUMBER,  WHILE  M  ALWAYS 

C  DENOTES  A  CELL  BOUNDARY  NUMBER.  SIMILARLY,  U,R,P,E  DENOTE  AVERAGE 
C  CONDITIONS  IN  A  CELL,  WHILE  CU,CR,CP,CE  DENOTE  CONDITIONS  AT  A  CELL 
C  BOUNDARY,  BEHIND  THE  WAVES  RESULTING  FROM  THE  DISCONTINUITY  THERE. 

C 

C  ********************************************************************* 

c 

c  INITIAL  AND  BOUNDARY  CONDITIONS 
C 

C  ********************************************************************** 
C 

C  INITIAL  CELL  PROPERTIES 
C 

C  DRIVER  SECTION  PROPERTIES 

DATA  (U(I,0),I-I,MD)/MD*0./ 

DATA  (R(I ,0) , 1-1 ,MD)/MD*5 . 66/ 


n  n  n  ono  o  no 


DATA  (P(I,0) ,I-1,MD)/MD*.487E6/ 

DATA  (E(I,0) , 1-1 ,MD)/MD* . 2151065E6/ 

DRIVEN  SECTION  PROPERTIES 
DATA  (U(I,0) , I-MD+1, MT)/MR*0./ 

DATA  (R(I , 0) , I-MD+1 ,MT) /MR* 1 . 19/ 

DATA  (P(I,0) , I-MD+1, MT)/MR*.102E6/ 
DATA  (E(I , 0) , I-MD+1 ,MT) /MR* . 2I51065E6/ 


NUM-0 
L-. FALSE. 

10  CONTINUE 
1-0 

WMAX-0. 

IF  (L)  NUM-NUM+1 
L-. FALSE. 

RNUM-NUM 
Z-RNUM/1000 . 

IF  (TIME(T) .GE.Z)  L- . TRUE . 

CONDITIONS  AT  THE  LEFT  BOUNDARY  (WALL)  ARE  REFLECTION  CONDITIONS. 

U(0,0)— -U(1,0) 

R(0 , 0)-R(l , 0) 

P(0,0)-P(1,0) 

E(0,0)-E(1,0) 

CONDITIONS  AT  THE  RIGHT  BOUNDARY  (WALL)  ARE  REFLECTION  CONDITIONS. 

U(MT+1,0)— U(MT,0) 

R(MT+1,0)-R(MT,0) 

P(,MT+i ,  0)-F(MT, U) 

E(MT+1,0)-E(MT,0) 

C 

IF  ((MOD(T,50).EQ.O).OR.(T.LE.5))  THEN 
WRITE  (5,100)  TIME(T) 

WRITE  (6.100)  TIME(T) 

CC  WRITE  (6,103)  T,CU(0,0)  ,T,CR(0,0)  ,T,CP(0  n^x.CFfO  O') 

WRITE  (5,102)  I ,T ,U(I , 0) , I ,T ,R(I , 0) ,I,T,P(I,0) , I ,T , E(I , 0) 

END  IF 
C 
C 
C 

C  **********************★************★***+***************************■*** 
C 

C  COMPUTE  THE  GENERAL  CONDITIONS  BEHIND  THE  WAVES. 

C 

Q  **************************■*******★****■**********************•**■*•*■*'***'*'* 

c 

C  GENERAL  NOTE:  DUE  TO  THE  APPROXIMATION  OF  A  FINITE  NUMBER  OF  CELLS 
C  WITH  CONSTANT  PROPERTIES  IN  EACH,  DISCONTINUITIES  EXIST  IN  THE 
C  PROPERTIES  AT  THE  CELL  BOUNDARIES.  THE  DISCONTINUITIES  ARE  RESOLVED 
C  BY  WAVES  (COMPRESSION  OR  EXPANSION)  THAT  ORIGINATE  AT  THE  BOUNDARIES. 
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o  ci  n  o 


C  THE  CONDITIONS  BEHIND  THESE  WAVES  ARE  THE  PRESSURE  (PO)  AND  VELOCITY 
C  (UO)  ON  THE  CONTACT  DISCONTINUITY,  AND  THE  DENSITY  ON  EITHER  SIDE  OF 
C  THE  CONTACT  DISCONTINUITY  (R2  AND  R3)  . 

C 

C 

C  COMPUTE  PO  AND  UO  BY  APPROXIMATE  FORMULAE  FOR  WEAK  WAVES. 

C 

1  IF  (ABS (P(I+1 ,0)-P(I,0)) .GT. (P(I+1 , U)/IOO . ) )  GO  TO  2 
IF  (ABS(U(I+1 ,0)-U(I,0)).GT.ABS(U(I+l,0)/100.))  GO  TO  2 
A-SQRT(0 . 25*G*(P(I , 0)+P( 1+1 ,0))*(R(I,0)+R( 1+1,0))) 

B— A 

PO— 0 . 5*(P(I ,0)+P(I+l , 0) )+0 . 5*A*(U(I ,0)-U(I+l,0)) 

UO-O . 5*(U(I,0)+U(I+1,0))+0.5*(P(I,0)-P(I+1,0))/A 
GO  TO  6 


COMPUTE  PO  AND  UO  BY  ITERATION  OF  FULL  EQUATIONS  FOR  STRONG  WAVES. 

2  PCD(0)-(P(I , 0)+P(I+l , 0) )/2 . 

N-0 

8  IF  (N.GT.50)  THEN 

PRINT*,  'ITERATIONS  FAILED  TO  CONVERGE  IN  50  TRIALS’ 

PRINT*,  '  N— ' , N , '  I-', I,'  T-’,T 

PRINT*,  ’  ' 

PO-PCD(N) 

GO  TO  9 
END  IF 

S— SQRi ( G*P ( I , 0 ) *R (1,0)) 

SI— SQRT(G*P(I+1 , 0)*R( 1+1 , 0) ) 

PI— PCD(N)/P(I ,0) 

PI1— PCD(N)/P(I+1 , 0) 

IF  ( (PCD(N)+I .) . GE . P( I , 0) )  THEN 
A— S*SQRT ( 1 . -G2*(l . -PI)) 

ELSE 

A— G1*S*(1  -PI)/(1. - (PI**G1) ) 

END  IF 

IF  ( (PCD(N)+1 . ) . GE . P(I+1 , 0) )  THEN 
B— S1*SQRT ( 1 . -G2*( 1 . -PT1) ) 

ELSE 

B— G1*S1*(1 . -PI1)/(1. - (PI1**G1) ) 

END  IF 

F— (B*P(I , 0)+A*P(I+l , 0)+A*B*(U(I , 0) -U(I+1 ,0) ) )/(A+B) 
FPCD(N)-MAX(PMIN, F) 

IF  (N.EQ.O)  THEN 
PCD ( N+l ) — FPCD ( N ) 

N-N+I 
GO  TO  8 
END  IF 

IF  ((FPCD(N-l) .EQ.PMIN) .AND. (FPCD(N) .EQ.PMIN))  THEN 

PO-PMIN 

PRINT*,  'NEGATIVE  PRESSURE  IN  ITERATION  WAS  RESET  TO  PMIN . ' 
PRINT*,  '  N-' , N , '  I-', I,'  T-',T 

PRINT*,  '  ' 

GO  TO  9 
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END  IF 

DPCD— PCD(N) -PCD(N-L) 

DFPCD-FPCD(N) -FPCD(N-l) 

IF  ( ABS (DPCD) . LT . ABS (PCD (N) /1000 . ) )  THEN 
PO-PCD(N) 

ELSE 

PCD(N+1)— (FPCD(N)*DPCD- PCD (N)*DFPCD)/( DPCD- DFPCD) 

N-N+l 
GO  TO  8 
END  IF 

9  U0-(A*U(I,0)+B*U(I+1,0)+P(I ,0)-P(I+I,0))/(A+B) 

C 

C 

C  COMPUTE  DENSITY  ON  EACH  SIDE  OF  THE  CONTACT  DISCONTINUITY . 

C 

6  R2-((G+1. )*PO+(G-l.)*P( 1+1,0) )*R(I+l,0)/( (G-l. )*P0+(G+1. )*P( 1+1,0) 
+  ) 

R3-((G+1.)*P0+(G-1.)*P(I,0))*R(I,0)/((G-1. )*P0+(C+1. )*P(I,0)) 

C 

C 

C 

C  ********************************************+*************■*''*'++*'**'**■**■* 

c 

C  COMPUTE  CONDITIONS  BEHIND  WAVES  AT  THE  LOCATION  OF  THE  CELL  BOUNDARY 
C 

C  *t+******************************************************************x-* 

c 

M-I 

WL—U( I , 0) - (A/R( I , 0) ) 

WR-U(I+1 ,  0)  +  (B/R(I+l , 0) ) 

WMAX-MAX (UMAX, ABS (WL) ,ABS(WR)) 

IF( ( (WL.GT.O. ) .AND. (WR . GT . 0 . ) ) . OR . ( (WL. LT . 0 . ) .AND. (WR.LT.O. ) ) )THEN 
IF  (WL.GT.O.)  THEN 
CU(M,0)-U(I,0) 

CR(M,0)— R(I ,0) 

CP(M, 0)-P(I ,0) 

ELSE 

CU(M,0)-U (1+1,0) 

CR (M , 0 ) -R ( 1+1 , 0 ) 

CP(M ,  0)-P( 1+1 ,0) 

END  IF 
ELSE 

IF  (UO.GT.O. )  THEN 
CU(M,0)-U0 
CR(M,0)-R3 
CP(M,0)-P0 
ELSE 

CU(M,0)-TJ0 
CR(M,0)-R2 
CP(M,0)-P0 
END  IF 
END  IF 

CE(M, 0)— CP(M, 0)/( (G-l . )*CR(M, 0) ) 

IF  ( ( (MOD(T , 50) .EQ.O) .OR. (T.LE. 3) .OR. (T.EQ. 399)) .AND. ( (MOD (M, 10) .E 
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+Q.0) .OR. (M.LT.10)))  THEN 

WRITE  (6,101)  M,T,CU(M,0) ,M,T,CR(M,0) ,M,T,CP(M,0) ,M,T,CE(M,0) 

END  IF 
I-I+l 

IF  (I.LE.MT)  GO  TO  1 
C 
C 
C 

Q  ************  1cic*ic*ic*ic**1c**ic*1cic**ic**ic*icicit*icicicicic*ic1c*icie'k'k*'k'k'k'k'k-k'k'k'k'k-k'k'k-k'k'k 

c 

C  COMPUTE  LENGTH  OF  CURRENT  TIME  STEP. 

C 

C  ********************************************************************* 
n 

IF (WMAX. EQ . 0 . )  PRINT*,  'WMAX-O.' 

TAU- .  9*H/VMAX 

TIME (T+l )— TIME (T)+TAU 

WRITE  (8,105)  T,TAU,TIME(T) 

C 

C 

C 

C  ********************************************************************* 

c 

C  UPDATE  CELL  PROPERTIES  USING  CONSERVATION  EQUATIONS. 

C 

C  ********************************************************************* 

c 

C  GENERAL  NOTE:  THE  FLUX  OF  PROPERTIES  INTO  EACH  CELL  IS  DETERMINED 
C  BY  CONDITIONS  AT  THE  CELL  BOUNDARIES  DURING  THE  TIME  STEP,  COMPUTED 
C  PREVIOUSLY. 

C 

1-1 
3  M-I 

R( I , I)— R( I , 0) - (TAU/H)*(CR(M , 0)*CU (M , 0) - CR(M- 1 , 0) *CU (M- 1,0)) 

MOM1-R (I,0)*U(I,0) 

MOM2-CP(M,0)+CR(M,0)*(CU(M,0)**2. ) 

MOM3-CP(M- 1 , 0)+CR(M-l , 0)*(CU(M- 1 . 0)**2 . ) 

U(I , 1)-(M0M1- (TAU/H)*(M0M2 -MOM3) )/R(I , 1) 

NRG1— R(M, 0)*(E(M,0)+(U(M, 0)**2 . )/2. ) 

NRG2-CR(M, 0)*CU(M, 0)*(CE(M, 0)+(CU (M , 0)**2 . )/2.+CP(M , 0)/CR(M , 0) ) 
NRG3— CR(M- 1 , 0)*CU (M- 1 , 0)*( (CE(M-1 , 0)+(CU(M- 1 , 0)**2 . )/2 . )+CP(M- 1,0) 
+/CR(M- 1,0)) 

E(I , l)-( (NRG1- (TAU/H)*(NRG2-NRG3) )/R(I,l))-(U(I,l)**2. )/2 . 
P(I,l)-(G-l.)*E(I,l)*Ra,l) 

IF  ( ( (MOD(T , 50) .EQ.O) .OR. (T . LE . 5) ) . AND . ((MOD(I , 10) . EQ . 0) . OR . ( I . LT . 
+10)))  THEN 

WRITE  (5,102)  I ,T ,U( I , 0) ,I,T,R(I,0) , I , T , P( I , 0) , I , T , E( I , 0) 

END  IF 
Y-I 

IF  (I.EQ.STA1)  THEN 

WRITE  (17,104)  Y, U(I , 0) , R( I , 0) , P( I , 0) , E( I , 0) , TIME(T) 

END  IF 
Y-I 

IF  (I.EQ.STA2)  THEN 


99 


nonn 


WRITE  (27,104)  Y,U(I , 0) , R( I , 0) , P(I , 0) , E(I , 0) , TIME(T) 

END  IF 

IF  ( (M0D(T, 75) . EQ . 0) .AND . ( (MOD (I ,10) . EQ . 0) .OR. (I.EQ.l)))  THEN 
IF((T.EQ.O) .OR. (T.EQ.55) .OR. (T.EQ.109) .OR. (T.EQ.164) .OR. (T.EQ.218) 
+.OR. (T.EQ.273) .OR. (T.EQ.328) . OR . (T . EQ. 382) . OR. (T.EQ.437) 

+.OR. (T.EQ.492) .OR. (T.EQ.546) .OR. (T. EQ . 598) . OR. (T.EQ.648))  THEN 
IF  (L)  THEN 

IF  (I.EQ.l)  WRITE  (7,106)  TIME(T) 

RI-I 

X-RI/100 . 

WRITE  (7,104)  X,U(I,0)  ,R(1 ,0) , P(1 ,0) , E(1 ,0) ,TIME(T) 

END  IF 
I-I+l 

IF  (I.LE.MT)  GO  TO  3 
C 
C 
C 

C  ********************************************************************** 

c 

C  UPDATE  CELL  PROPERTIES  AND  PROCEED  TO  NEXT  TIME  STEP. 

C 

C  ********************************************************************** 

c 

DO  7,  I-l.MT 
u ( I , 0 ) —U (1,1) 

R(I,0)-R(I,1) 

P(I,0)-P(I,1) 

E(I,0)-E(I,1) 

7  CONTINUE 
T-T+l 

IF  (T.LE.TT)  GO  TO  10 
C 
C 
C 

C  ********************************************************************** 

c 

C  FORMAT  STATEMENTS 
C 

C  ********************************************************************** 

c 

100  FORMAT  (1H1, 'TIME-' , F8 . 6/) 

101  FORMAT ( 5X , 'CU( ' , 13 , ' , ' , 13 , ' )-' , E16 . 9 , 6X, ' CR( ',13 , ' , ' , 13 , ' )-' , E16 . 9 
+,6X,'CP(',I3,',',I3,')-',E16.9,6X,’CE(',I3,',',I3,')=-', E16 . 9) 

102  FORMAT ( 5X , 'U( ',13 , ' , ',13 , ' )-' , F16 . 9 , 7X, ' R( ’ , 13 , ’ , ',13 , ' )-' , E16 . 9 , 
+7X, ' P( ' , 13 , ' , ' ,13, ')— ' , E16 . 9 , 7X, 'E(' ,13, ' , ' ,13, ')—' , E16 . 9) 

103  FORMAT  (5X,'CU(  0 , ' , 13 , ' )-' , E16 . 9 , 6X, ' CR(  0 , ',13 , ' )-' , E16 . 9 , 

+6X, ' CP(  0, ' ,13, ')-' ,E16.9,6X, 'CE(  0 , ' , 13 , ' )-' , E16 . 9) 

~U4  FORMAT  (IX, 6(E14. 7) ) 

105  FORMAT  (IX, 110 , 2F16 . 10) 

106  FORMAT  (1X.E14.7) 

END 
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